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Chapter 1 
Introduction 



1.0.1 Motivations 

Strongly interacting quantum many-body systems have been one of the main chal- 
lenges of quantum physics and are still not well understood in many aspects; many 
novel intriguing phenomena may in fact be originated from the strong interactions 
among particles in these systems [T]. A strongly interacting system can be described as 
a system where one can not define a small parameter on which a perturbative theory 
can be built. This complication inspired the development of numerical approaches 
based on the variational principle [2] and also quantum simulations [3l H] that, in the 
case of bosonic systems, are in principle "exact" . In this work we have considered 
two systems that can be regarded as the archetype for neutral strongly interacting 
systems: 4 He, which is a bosonic system, and its fermionic counterpart, 3 He. More 
specifically, in this work we employed Quantum Monte Carlo (QMC) techniques at 
zero and at finite temperature, respectively the Path Integral Ground State [3] (PIGS) 
and the Path Integral Monte Carlo[4J (PIMC), to the study of a system of two di- 
mensional 3 He (2<i- 3 He) and to the study of 4 He adsorbed on Graphene-Fluoride (GF, 
called also Fluorographene) and Graphane (GH), namely two corrugated substrates 
that can be derived[6] from Graphene. Our main purpose in the case of adsorbed 4 He 
was the research of new physical phenomena, whereas in the case of 2d 3 He it was 
the application of novel methodologies [7] for the study of static properties of Fermi 
systems and the extension of such methodologies for an ab-initio study of the low 
energy excitations of a strongly interacting fermionic system. 

Apart from being both strongly interacting, the systems that we have considered 
are interesting also from a methodological point of view, as they can be used to test the 
limits of the employed techniques. In the case of 2d- 3 He the main technical difficulty 
relies in the well known sign problem^, which, on one side, poses a severe limit on 
the number of particles that can be simulated by QMC and, on the other side, limits 
the study of imaginary-time dynamics to small values of imaginary-time. For 4 He on 
GF and GH the geometry of the confinement gives rise to rare tunneling events that 
are relevant in both the static and dynamic properties of the system and must thus 
be correctly described by the used QMC technique. The relevance of these system is 
also increased by the fact that experiments are feasible on both systems, indeed for 
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2d- 3 He there is already a number of experimental works in literature and comparison 
with experiments has been done wherever possible; our study of 4 He on GF and GH 
instead is novel and up to now there is no experimental data with which we can 
compare our predictions, however it has been shown in Ref. [5] that the substrates 
that were considered are available to the experimentalists and we hope that this work 
will inspire some new experiments on this topic. 

In the remnant of this section we give a first introduction of the systems that have 
been considered. 

Two dimensional 3 He 

Two dimensional bulk 3 He at zero temperature is a model well suited for the study of 
strongly correlated Fermi systems. This is because, as shown in Ref. [E], the model is a 
good approximation for liquid 3 He adsorbed on preplated graphite substrates. Indeed, 
much experimental work has been done on such systems, we mention heat capacity 
measurements in Ref. [HJ [TU] and more recently [TT] , the study of the thermodynamic 
behavior of the second layer of 3 He has been done in Ref. |12j . the study of magnetic 
properties of liquid 3 He films [T31 [H] and the study of low energy excitations with 
neutron scattering experiments [15] IT6] ; another feature of such systems is also the 
possibility to realize small clusters with a controlled number of particles [17]; this 
is appealing because those systems can possibly be simulated with "exact" QMC 
techniques. 

Also from the theoretical side, 2d 3 He has been the subject of many works, we 
mention the thermodynamic study of 2d Fermi liquid with and without external 
magnetic field [T8l IT9] . a many-body study of elementary excitations is reported in 
Ref. [20J, a QMC computation of the zero temperature equation state of pure 2d 
3 He[2T] and an estimation of its effective mass[22j; 

The experimental works in Ref. [8] revealed that quasi-two-dimensional 3 He has 
a nearly perfect Fermi liquid behavior, in particular, they showed that the effective 
mass m* and the spin susceptibility x/Xo increase with the density. This behavior, 
consistent with a divergence of m* near the freezing density, has been interpreted [23] 
as a signal of Mott transition to an insulating crystal. On the other hand, quasi-two- 
dimensional 3 He has been studied by theoretical means [22] that suggested that the 
freezing and the divergence of m* may not have the same physical origin, in particular 
the freezing density is influenced by the preplated substrate. In this context, the study 
of the strictly 2d 3 He becomes valuable in order to isolate the effect of correlations 
on the system near freezing density. A further advantage in the theoretical study of 
this system is that the properties of the liquid phase are largely independent on the 
choice of the substrate and thus it is possible to make a comparison with experimental 
data[SJ. An even greater interest in 2d 3 He has been also inspired by the recent 
work in Ref. [201 121] in which, for the first time, the collective zero-sound mode has 
been observed as a well defined excitation crossing and possibly reemerging from the 
particle-hole continuum. 

We have thus performed a Quantum Monte Carlo study of a two-dimensional 
bulk sample of 3 He using the unbiased Fermionic Correlations (FC) technique that 
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has been successfully employed in the 2d electron gas in Ref. [TJ. This technique 
is a formally exact method that makes use of bosonic imaginary-time correlation 
functions of operators suitably chosen in order to extract fermionic energies. In this 
work we computed the energy per particle as function of the polarization of the system 
at different fluid densities, from this data we obtained a spin susceptibility that is in 
very good agreement with experiments. As a further study of the system, we have 
extended the FC method to study dynamical properties; we computed an ab-initio 
low-energy excitation spectrum of 2d 3 He obtaining a well defined zero-sound mode 
in remarkably good agreement with Ref. [24J. 

4 He on Graphane and Graphene-Fluoride 

Experiments on the adsorption of Helium on Graphite have been carried out in the 
seventies at the University of Washington; those experiments revealed for the first 
time a behavior corresponding to a two-dimensional gas. Moreover, the appearance 
of a peak in the specific heat of 4 He near a critical temperature T c = 3 K showed 
evidence of a phase transition from a high T fluid to a low T commensurate (\/3 x 
V3 R30°) phase, an ordered phase in which the 4 He atoms are localized on second- 
nearest neighbors hexagons. Following this discovery, a number of experimental and 
theoretical works followed and now the Helium monolayer on Graphite is probably 
one of the most studied adsorbed quantum systems. 

On the experimental side we mention specific heat measurements in Ref. [251 I2B] , 
chemical potential measurements in Ref. [27] and neutron scattering experiments 
in Ref. |28j. The phase diagram of the first layer of 4 He on Graphite has been 
inspected in Ref. [291 USB ED] . As for the second layer, we mention the experimental 
work in Ref. [32] . Superfluid properties of Helium on Graphite were investigated in 
Ref. [331 El E5]. 

On the theoretical side we mention the work on the interaction potential of He 
on Graphite by Carlos and Cole [36] and a study on the possible commensurate 
solid phases of the second layer presented in Ref. [371 138] . There are also many 
simulations [391 HQl SB S21 H3J E] on strictly 2d 4 He. As for Helium on Graphite, 
the role of corrugation has been studied with Path Integral in Ref. [13] whereas the 
properties of the adsorbed layers have been studied with Monte Carlo simulations in 
Ref. [HI HH1 H7_l HE] and more recently in Ref. [49] . There has also been works on He- 
lium on Graphene, the phase diagram has been calculated in Ref. [50] and superfluid 
properties in Ref. [5T] . 

The availability of Graphene and especially its derivatives like Graphane and 
Graphene-Fluoride makes possible the study of new adsorbed systems. No special 
phenomenon is expected for Helium adsorbed on Graphene because the interaction 
is geometrically similar to that on graphite, but in the case of GF and GH the ad- 
sorption potential is qualitatively different from the case of Graphite and indeed we 
found a unique behavior of the adsorption system. The difference of GF and GH from 
Graphite is due to their conformation; GF and GH are respectively Graphene sheets 
to which are chemically bonded planes of either Fluorine or Hydrogen atoms; in the 
case of GH, for example, the substrate is made of a Graphene sheet with Hydrogen 
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atoms attached above and below the C atoms, in an alternating pattern. Such atomic 
structure provides an Helium-substrate interaction potential which, compared with 
the Helium-Graphite potential, has twice the number of adsorption minima located 
on an honeycomb lattice; compared with Graphite, the tunneling between the ad- 
sorption sites of GF and GH is also enhanced along three spatial directions that cross 
saddle points of the potentials. These properties of the GF(GH) adsorption poten- 
tial, as shown in Sec. |5j not only confine Helium in a multi-connected space but also 
destabilize the analogue of the v3 x a/3 R30° on Graphite: we found that the ground 
state at equilibrium density, for both GF and GH, is indeed a modulated superfluid 
that in GF has anisotropic rotons in the excitation spectrum. Also high coverages 
of 4 He monolayer on GF and GH show novel properties that have been described 
in Sec. [5] we found in fact a stable commensurate solid phase that is the analogue 
of the theoretically predicted 4/7 phase on Graphite, moreover we have preliminary 
evidence that this solid phase possesses also a relevant superfluid fraction. 

1.1 Implemented Methodologies 

Quantum Monte Carlo methods are largely employed in the study of strongly in- 
teracting quantum systems; the main reason for that is because they can provide 
expectation values that can be in principle "exact" in the case of Bose systems. In 
the case of Fermi systems, QMC methods are still an highly accurate tool. The word 
"exact" here means that the used approximations may be reduced below the statis- 
tical error of the QMC method. To make a few examples of successful applications 
of QMC methods, we mention the quantitative evaluation [52] of the Bose-Einstein 
condensate fraction in liquid 4 He at zero temperature, the phase diagram of 4 He ad- 
sorbed on Graphite [49J and, more recently, the low energy excitation spectrum [53] 
of 4 He at zero temperature and the computation of the normal-state equation of a 
Fermi ultra-cold gas at unitary regime [54] . 

The first QMC method that appeared was a variational technique named Vari- 
ational Monte Carlo[2] (VMC). This technique expresses a zero temperature expec- 
tation value on a given family of variational wave functions as a multi-dimensional 
integral and then compute the integral with the Metropolis algorithm [55]. Originally 
it was implemented with Jastrow wave functions [56], but better classes of trial wave 
functions were introduced; it is worth to mention here the Shadow Wave Functions 
(SWF) for Bosons[57] and for Fermions[58j. that introduce many-bodies correlations 
in an implicit way and is able to describe a system in both the liquid and the solid 
phases, without introducing explicitly any equilibrium lattice for the solid state. Be- 
yond the variational level, the first introduced "exact" QMC technique was the Dif- 
fusion Monte Carlo [59J (DMC) that solves the Schrodinger equation for the ground 
state of a many-body system taking advantage of its similarity with the diffusion 
equation in imaginary time. Another exact technique valid at zero temperature that 
was developed soon after DMC is the Green's Function Monte Carlo [60] (GFMC); 
this method exploits an integral formulation of the Schrodinger equation in order to 
express ground state quantum averages; on the same line, another very successful 
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method is the Path Integral Ground State [3] (PIGS) that expresses a ground state 
expectation value through Feynman's path integrals as a sufficiently long imaginary- 
time evolution of a trial wave function; an improvement of DMC that had been 
introduced in the same years of PIGS is the Reptation Monte Carlo [6"T] (RMC). Like 
in the case of VMC, better trial wave functions have been constantly introduced in 
PIGS; one of the last advancements in zero temperature path integral simulations on 
Bose systems is the Shadow Path Integral Ground State [62J (SPIGS) which makes 
use of SWF as trial wave functions. A very strong feature of PIGS and SPIGS is 
that they are formally similar to the Path Integral Monte Carlo[4j (PIMC) method; 
PIMC, in fact, uses Feynman's path integrals in order to compute quantum thermal 
averages; apart from that, its remarkable formal similarity with PIGS comes also from 
the similarity between the thermal density matrix and the quantum imaginary-time 
evolution operator. This feature has a practical value because the two methodologies 
can be implemented within the same framework. 

The mentioned "exact" methodologies, if applied to Fermi systems, suffer from 
the sign problem^. This problem occurs because the Fermi symmetry introduces a 
nodal surface in the ground state wave function (or in the density matrix elements in 
the case of PIMC) that, as consequence, is no longer a probability density that can 
be sampled with Monte Carlo; the same problem is also present in Bosonic systems if 
an excited state instead of the ground state is considered. There are workarounds but 
they result in a signal to noise ratio that decreases exponentially with the number of 
particles; exact Fermi simulations, as well as the study of the excitations of Bosonic 
systems, are thus restricted to system with small number of particles. Among the 
adaptations that allow the QMC computation on Fermi systems there is the Fixed 
Node[63j approximation (FN), a variational technique that approximates the true 
nodal surface of the ground state with that of a trial wave function, and its evolution, 
the Released Node [61], that has shown to be exact for small systems [61] : we also 
mention the Restricted Path[65] method that extends PIMC to Fermi systems and 
a more recent evolution [66] of the DMC method that gives exact results for small 
systems. 

The mentioned techniques work in real coordinates space; another rather new 
and promising approach to the study of Fermi systems is the formulation of novel 
QMC techniques; we mention here the Auxiliary Fields Quantum Monte Carlo [67] 
(AFQMC) and the Bold Diagrammatic Monte Carlo [SH] (BDMC). 

In this work we have studied Fermi systems with SPIGS; for this purpose we 
have adopted another recently developed technique named Fermionic Correlations [7] 
(FC). FC can be defined as a "cross-over" technique because its basic idea is to 
obtain informations on a Fermi system through the computation of an imaginary- 
time correlation function on a fictitious Bose system; with this approach, the sign 
problem is avoided and the simulation is in principle exact. However, to obtain the 
informations on the Fermi system from the imaginary-time correlation function one 
has to compute a numerical inversion of the Laplace transform in ill-posed conditions, 
this is a difficult inverse problem that, again, results in severe limits on the number 
of particles that can be studied. If the number of particles is small enough, however, 
the FC technique is an unbiased, ab-initio method that gives access to the energy 
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(and possibly its derivatives) of strong interacting Fermi systems; moreover, the FC 
technique has been extended in this work to study collective excitations of Fermi 
systems. 

1.2 Thesis Outline 

In this work we have made an "unconventional" choice: instead of making a single 
chapter devoted to the full theoretical introduction of the methodologies, we have 
introduced the essentials in chapter [2] and the technical details of the methodologies 
in the chapter after the conclusions. With this choice, a reader that is not interested 
in technical details can safely ignore anything written after the conclusions. 
This document is organized as follows. 

• The present section provides a background on both the studied physical systems 
and the employed methodologies. 

• In chapter [2] we provide a basic description of the PIGS and PIMC techniques. 
In this chapter, a methodological work is also presented. We show that, on a 
realistic model potential for 4 He, the PIGS method does not suffer from any 
bias deriving from the choice of the trial wave function. 

This work has been published on J. Chem. Phys., 131, 154108 (2009). 

• In chapter [3] is presented the study of 2d 3 He at zero temperature with the FC 
technique. The energy of the system for various densities and polarizations is 
reported as well as the resulting spin susceptibility as function of the density. 

The work includes also comparison with experimental data and Fixed Node 
simulations and has been published on Phys. Rev. B, 85, 184401 (2012). 

• In chapter [I] we adapted the FC technique to study the excitations of a Fermi 
system. The reader can find an ab-initio computation of the dynamic structure 
factor of 2d 3 He at zero temperature compared with recent experimental data, 
the static response function and the approximate static structure factor. 

These results are in preparation for submission to Phys. Rev. B. 

• In chapter [5] we present the study of Helium adsorbed on Graphene-Fluoride 
(GF) and Graphane (GH). The section will present one body properties, such 
as the ground state energy of one atom of 3 He and 4 He on GF and GH and 
the first energy band in the four cases; it will treat then many-body properties 
of the first layer of 4 He, such as the stability of various commensurate phases, 
the equation of state at zero temperature, the condensate fraction in the liquid 
phases, the zero temperature low energy excitation spectrum at the equilibrium 
density and superfluid properties at both zero and finite temperature. We also 
present preliminary data on a possible supersolid phase present at high coverages 
on both GF and GH. 
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Many of these results have been published on: 

J. Phys.: Conference Series 400, 012010 (2012) - proceedings of the LT26 
conference. 

J. Low. Temp. Phys. - proceedings of the QFS2012 conference. DOI: 
10.1007/sl0909-012-0770-9 

Phys. Rev. B. 86, 174509 (2012). 

• In chapter [6] we draw the conclusions of this work. 

• in chapter [7] the computational details of the PIGS and PIMC methods are 
thoroughly described, from the mathematics of the Markov chain to the im- 
plementation of the Metropolis algorithm and the derivation of estimators that 
compute expectation values of various physical quantities. 
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Chapter 2 

Path Integral Methods 



In this Chapter the general basis of two Monte Carlo techniques will be described; 
the technical details will instead be discussed in Chapter [7} The method used for zero 
temperature simulations is the Path Integral Ground State [1] (PIGS) whereas, that 
used for finite temperature simulations is the Path Integral Monte Carlo [2] (PIMC). 
The PIGS and the PIMC techniques are "exact" methods if the studied system has 
the Bose symmetry; the word "exact" in the context of Quantum Monte Carlo (QMC) 
means that the systematic errors due to the used approximations can be arbitrarily 
reduced below the Monte Carlo statistical uncertainty The two techniques have 
also a similar formalism. For this reason, they are easily implementable in a unified 
computer library. 



2.1 Path Integral Ground State 

In Sec. ITlwe show that, using Monte Carlo techniques, it is indeed possible to sample 
an arbitrary probability distribution and that with the resulting sampling it is possible 
to evaluate ^-dimensional integrals. We now specialize that methodology to the 
problem of calculating the expectation values of a bosonic iV-particle system. 

Let's thus consider a system of iV atoms of mass m at a temperature T = K, in 
a box of volume V& in periodic boundaries conditions, with an interatomic potential 
V(r), the Hamiltonian operator is 

H = f + V (2.1) 

where the kinetic term is 

fc2 N 
1=1 

and the potential term is 

V = J2v(\f t -f 3 \) (2.3) 

i<j 
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this Hamiltonian is used for ease of writing, but a more general Hamiltonian with 
anisotropic interactions and external potential can be used as well. Defined ^{R) as 
the ground state wave function, we want to compute the quantity 

(6) = J dRO(R)^ 2 (R) (2.4) 



where R = {r{\^ =l is a many-body variable and r*j is the position of the i-th particle 
of the system and O is an operator that is diagonal in the coordinate representation. 
The square of the wave function, \1/ 2 (_R), real and nodeless because we are considering 
a bosonic system, is proportional to the the probability distribution to be sampled 
with the Metropolis algorithm (see Chapter [7]). The quantity \l/ 2 (i?) is in general 
unknown but a workaround that has been very successful among T=0 K methods is 
to exploit the quantum evolution in imaginary time. 

Given an initial state (0)), the quantum time-evolution is determined by the 
Schrodinger's equation and 

|tf (*)) = e -K*# |tf (0)) (2.5) 

where the time evolution operator U(t) = e~^ tH . If is an eigenvector of H, its 
overlap with the state \^(r)) can be expressed as 



(*,|* (r)> = {^\e~ rH \^) (0)) (2.6) 

3 

where we have defined the quantum imaginary-time evolution operator U (r) = e~ rH 



by substituting r = \t. Eq. ( |2~1?] ) can be rewritten as ($j|\l>(r)) = e~ rEl ($»|* (0)). 
For a sufficiently long r, if the initial state ^ (0) is not orthogonal to the ground state, 
only the eigenstate corresponding to the lowest eigenvalue has a relevant overlap 
on the evolved trial wave function |\l/(r)). The ground state wave function \l/o m 
coordinate representation can be thus expressed as the r — > oo limit of an imaginary 
time evolution of an arbitrary trial wave function v&r provided that (^qI^t) 7^ 



e -T{H-E 0) y 



% = lim - - ; T . (2.7) 

The normalization factor is not involved in the Monte Carlo sampling; within the 
Green's function formalism, the ground state wave function can be approximated 
with ^r(R), 



q T (R) = _1 J dR'G {R, R', r) ^ T (R) (2.8) 

where M is the normalization constant and the term G(R, R ; , r) = (R\e~ rH \R') is 
the Green's function or density matrix. Here, the expectation value has just merely 
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been rewritten in term of the Green's function, but the Green's function for a suf- 
ficiently large r is still a generally unknown quantity. There are, however, known 
analytic approximations of the Green's function that are valid for small imaginary 
time 6~t and the Path Integral formalism provides a way to express a large r Green's 
function as a convolution of smaller imaginary-time Green's functions. This comes 
from an important property of the density matrix, 



-tH 



-StH 



m 



(2.9) 



where 6~t 



M ' 



In the coordinate representation, the product becomes a convolution 



» . M-l 

G(R 1 ,R m+u t) = J ... J dR 2 ...dR M Y[G(R v R J+1 ,5r 



(2.10) 



A density matrix at imaginary time r can be represented as a convolution of M density 
matrices at smaller imaginary time t/M. This convolution is the Path Integral and, 
as the name Path Integral Ground State may suggest, it is a fundamental element for 
the quantum simulation techniques that have been used throughout this work. 



Combining Eq. (2.4) with (2.7) and (2.8), a quantum average on the ground state 
thus becomes 



d ) = Jff (j\ dR ^j ^t{Ri)0{R m/2 ) 

M-l 

x YlGiR^R^d^^TiRM) 



X 



(2-11) 



In the case that r/2 is sufficiently large to have a good approximation of Eq. (2.7), if 
the operator O commutes with e~ rH , then it can be applied at any position k of the 
path integral and it will give the ground state expectation value, if otherwise, [O, H] ^ 
0, then the operator applied at positions k = 1 and k = M will give mixed expectation 
values (* T |O|* ) and for k = 2...M/2 the expectation values (^ T |e- (fc - 1)<5T ^|O|*(0)) 
will converge to the ground state value. To go further and obtain an explicit definition 
of the quantum expectation value, an analytic approximation of the small imaginary- 
time Green's function must be used. The simplest one is the Primitive Approximation 
(PA); more advanced approximations are illustrated in Appendix. [Bj The PA consists 
in neglecting the commutator between T and V when factorizing the density matrix 
e~ 5rH , the error associated with this approximation is of the order St 2 . 



e SrH _ e SrT e - 



StV 



(2.12) 
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In this approximation the matrix elements of the two factors are easy to obtain, 



R, 



-StT 



R; 



Ri+1 
8tV 



1 



4At 



(AnXdr) 2 
R i+1 ) = e- 5TV ^6(R l ,R i+1 ) 



where we have defined A 



2m' 



\Ra — R 



»+i| 



2^7 = 1 



i+1 



(2.13) 
(2.14) 

1 2 

and V (R m ) = 



I i j 



Si<i v ( r «7) ' ^ e use ^ e snor thand notation r-j 

The ground state expectation value of an operator O that is diagonal in the 
coordinate representation becomes 



O 



1 

Jf 



M-l 



\ [ dRi * T (Ri)e 



-V(Ri) f 



e -^v(R i+1 ) ^ Rm/2) ^(^12.15) 



i=l 



where the primitive approximation has been written in a symmetric form, namely 



G PA (Ri,R i+1 ,ST) = e 



-V(Ri) f 



{Rj-R i+1 y 



e -4fV(fli+i) 



In the limit M — > oo, Eq. (2.15) becomes exact due to the Trotter formula : 

M 



-tH 



lim 

M— yoo 



,(StT) J-StV) 



(2.16) 



(2.17) 



Chosen a sufficiently high M, then, the error in Eq. (2.15) can be arbitrarily reduced, 
and, chosen a sufficiently high r, an arbitrary precise description of the ground state 
can be obtained. The integral can be evaluated with Monte Carlo and the multi- 
dimensional probability distribution to sample with the Metropolis algorithm is 



A/-1 



P ({R n }) = *T(Rl) II G ^ G ^' *t(Rm) 



(2.18) 



The value of r that is sufficiently high to have convergence depends by a good 
degree on the choice of the trial wave function ty?- A trial wave function that has an 
high overlap on the ground state could, in fact, enhance the convergence of Eq. (2.7). 
We emphasize that this, however, is not necessary to obtain unbiased results: a very 
strong feature of PIGS, as we have shown in Ref. [3], is that, indeed, the results of a 
PIGS calculation do not depend on the choice of the trial wave function. This is what 
we are going to show in Sec. \2.1.1 the practical implementation of the Metropolis 



algorithm will be instead described in Sec. 7.1.4 



2.1.1 Quantum— Classical Isomorphism 

Although path integrals and quantum evolution in imaginary time are very abstract 
topics, there is a simple interpretation of the probability distribution (2.18) that allows 
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for an easy visualization of the Metropolis sampling; besides the practical advantages, 
it is also an interesting example in which many aspects of both the mathematics of 
Markov chains and the physics of the system take life in a fictitious classical system 
made of beads that have very special interactions between each other; this is why this 
correspondence is called quantum-classical isomorphism. More specifically, Eq. (2.18) 
is the partition function of a classical system of N polymers composed by M beads 
that have special interactions. The kinetic term of the Hamiltonian represents the 
interaction between adjacent beads of the same polymer whereas the potential term 
of the Hamiltonian maps onto the interactions between beads of different polymers. 
A polymer is essentially a set of beads corresponding to some integration variables 
in Eq. (2.18), namely the z-th polymer is {f^ 5 } _ x . The length of the polymer is 

r; due to the analogy with the quantum evolution operator e~ %rB , this length is an 
imaginary time. The index j represents the position of the bead in the polymer, this 
position corresponds to a discrete imaginary time Tj = j5r = jjj- The discretization 
of the polymer in the imaginary time is named time-step. The mean square displace- 
ment of the beads in a polymer represents the indetermination of the position of the 
corresponding particle in the quantum system. A configuration of these special inter- 
acting polymers is thus defined by a set of coordinates {^*/}, where i represents the 
polymer and j represents the bead in the i-th polymer; Figure |2.1 shows a schematic 
representation of the polymers in PIGS and their correlations. A quantum observable 
is mapped to an operator that acts on such configurations that in this context shall 
be referred as estimator. 



Same polymer 
> = Bead 

O — I — O = Kinetic correlation 



= Potential correlation 



Same imaginary -time 

Figure 2.1: Schematic representation of the polymers in PIGS and their correlations. 



Here on, we will focus on Bose systems. In Sec. 7.1.4 the practical implementation 



of the Metropolis algorithm will be described; in that section, two different Monte 
Carlo algorithms used to sample the space of permutations will also be shown. 



Path Integral Ground State in action 

In my work for the Master degree I developed a library that can run PIGS simula- 
tions as well as Path Integral simulations at finite temperature (Path Integral Monte 



Carlo, PIMC, see Sec. 2.2). The work presented in this section is one of the early 
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employments of the library and it is a benchmark for both the library itself and the 
PIGS technique. A benchmark for the library because the library came through ex- 
tensive testing during this work, a benchmark for the PIGS technique because this 
work shows that PIGS is really unbiased, in other words, the choice of the trial wave 
function does not affect the final results, provided that the projection time is large 
enough and the time-step 5t is sufficiently small. This is a very strong feature of 
PIGS and is a necessary condition for a truly ab-initio method because it allows the 
study of a quantum system even if we don't know anything about its many-body 
ground state wave function. 



Test systems We have considered two bulk phases of a many-body strongly inter- 
acting Boson system: liquid and solid 4 He. Dealing with low temperature properties, 
4 He atoms are described as structureless zero-spin bosons, interacting through a re- 
alistic two-body potential, that we assume to be the HFDHE2 Aziz potential [3]; we 
remark here that our results are thus valid on this interaction potential and have not 
general validity. 

For the liquid phase, we have considered a cubic box with periodic boundary 
conditions, containing N = 64 atoms at the equilibrium density p\ = 0.0218A -3 . For 
the solid phase we have considered a cubic box with periodic boundary conditions 
designed to house a fee crystal of iV = 32 atoms at the density p s = 0.0313A -3 . In 
both cases we add standard tail corrections to the potential energy to account for the 
finite size of the system by assuming the medium homogeneous (i.e. g(r) = 1) beyond 
L/2, where L is the size of the box. Obviously, this is not an accurate assumption 
specially for the solid phase in such a small box, but our main purpose here is to show 
that PIGS method is able to reach the same results independently on the considered 
initial wave function. Computations of ground state properties of bulk 4 He with 
accurate tail corrections can be found in the current literature. [2J E] 



Trial wave functions The trial wave functions commonly used within the PIGS 
methodpQ are the variational Jastrow wave function (JWF) for the liquid and the 
Jastrow-Nosanow (J-NWF) for the solid. A JWF represents the simplest possible 
choice of wave function for strongly interacting Bosons [7] and it contains only two- 
body correlations. Using a McMillan pseudopotentialjS], the unnormalized JWF reads 
as 

N -ifb\ m 

^wfor)= n e 2{rijJ ■ ( 2 - 19 ) 
i<j=i 

The physical meaning of this JWF is that, due to the sharp repulsive part of the 
interaction potential V in the Hamiltonian H, 4 He atoms prefer to avoid each other. 



In the J-NWF the JWF is multiplied by a term like the one in Eq. (2.22) below, 
that localizes the particles in a crystalline order. In this work, however, in order 
to explore the convergence properties of the PIGS method, we have considered two 
wave functions of "opposite" quality: the best available one, that is the shadow wave 
function, and the poorest imaginable one, i.e. the constant wave function. As we 
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shall see, JWF will be considered only when computing the one-body density matrix 
in the liquid phase. 

The constant wave function is the ground state wave function of the ideal Bose 
gas, 

*Pcwf{R) = 1. (2.20) 

It carries no correlation at all. We choose this wave function because, allowing an 
unrestricted sampling of the full configurational space, it results in no importance 
sampling. Then the whole imaginary time projection procedure is driven only by the 
short imaginary time Green's function G(R, R' , St), without any input, and then any 
bias, from the initial state. Thus at the starting point the system is made up by free 
particles; if after a long enough imaginary time projection, PIGS turns out to be able 
to reach a strong correlated quantum liquid and quantum crystal by itself we can 
safely believe that no variational bias affects PIGS results. 

On the other hand, we choose as ipx a SWF optimized with a variational compu- 
tation in order to have as reference results the ones coming from the projection of an 
initial wave function that is more accurate as possible, i.e. from a wave function whose 
overlap with the exact ground state is known to be large. In the SWF, additional 
correlations besides the standard two body terms are introduced via auxiliary vari- 
ables which are integrated out[U]. This is done so efficiently that the crystalline phase 
emerges as a spontaneously broken symmetry process, induced by the inter-particles 
correlations as the density is increased, without the need of any a priori knowledge of 
the equilibrium positions and without losing the translationally invariant form of the 
wave function. Thus SWF is able to describe both the liquid and the solid phase with 
the same functional form and it is explicitly Bose symmetric. The standard SWF 
functional form reads 



VswfOR) = MR) J dSK(R,S)(j) s (S) (2.21) 
where S = (sj, s 2 , . . . , sjv) is the set of auxiliary shadow variables, <p r (R) is the stan- 



dard Jastrow two body correlation term (2.19), K(R,S) is a kernel coupling each 



shadow to the corresponding real variable, and ip s (S) is another Jastrow term describ- 
ing the inter-shadow correlations. Due to its analytical expression, the introduction 



of the SWF defined by Eq. (2.21 ) in a PIGS simulation consists in adding a timeslice 
at each extremity of every polymer. These newly added timeslices have special cor- 
relations; namely there are real-shadow intrapolymer correlations defined by K(R, S) 
and shadow-shadow interpolymer correlations defined by <j) s (S). As consequence, the 
PIGS method has to be extended with Metropolis moves that accordingly involve the 
introduced shadow timeslices. 

As usual [10], we take K(R, S) Gaussian and, as pseudopotential in <j) s (S), we use 
the He-He potential V rescaled in both amplitude and distances. The variational 
parameters we use were chosen in order to minimize the expectation value of the 
Hamiltonian H and are reported in Ref. |10j . Nowadays the SWF represents the 
best available variational wave function for 4 He systems. [6] Recently, it has been 
estimated [TT] that, when describing a two dimensional solid, the overlap of the SWF 
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with the true ground state is of about (0.998)^, which ensures a fast convergence rate 
when projected within the PIGS method. The properties of the SWF are so peculiar 
that the PIGS method that has a SWF as ipT deserves an its own name and is dubbed 
SPIGS: Shadow Path Integral Ground State method. [T2"| IT5] 

In order to test how robust PIGS is, we consider also a wave function that describes 
the wrong phase: for the liquid phase we consider a Gaussian wave function, where 
each particle is harmonically localized around fixed positions {roi} 

N 

^GWF(i?)=n e " cK " r " 12 ' ( 2 - 22 ) 

8=1 

i.e. ipT it the wave function of an Einstein harmonic solid. The parameter C = 
8 A~ 2 is arbitrary and it is was chosen to ensure a strong localization of the particles 
around the positions {rQi} that were taken over a regular cubic lattice within the 
simulation box. This wave function is evidently not translationally invariant and 
not Bose symmetric. Furthermore it does not contain any correlation between the 
particles, and all the information that it carries is that of a crystalline system, i.e. 
GWF is an extremely poor wave function for the liquid phase. This "bad" initial 
wave function will provide a stringent test on the convergence properties of the PIGS 
methods. 

As far as the one-body density matrix computation in the liquid phase is con- 
cerned, the values of the parameters b and m in the JWF have been chosen equal to 
the ones of the corresponding Jastrow term in the SWF. 



Small time Green's function One of the fundamental elements of path integral 
projection Monte Carlo methods is the imaginary time Green's function G(R, R', r), 
whose accuracy turns out to be crucial to the convergence to the exact results. The 
functional form of G for a generic r is unfortunately not known with exception of few 
particular cases, such as, for example, the free particle and the harmonic oscillator, 
but accurate approximations of G are obtainable in the small r regime j2j [HI [15]. In 
this work, we have chosen the Pair-Suzuki approximation[16j for the imaginary time 
propagator, which is a pair-approximation of the fourth-order Suzuki-Chin density 
matrix. [TJ] 

The Suzuki-Chin approximation is based on the following factorization of the 
density matrix: 

e -25rH „ e -fv ee -6rf e -^V Ce -6rf e -^ ^33) 

where T is the kinetic operator and V e and V c are given by 

X 2\ N 

K = ^+^B F i) 2 (2-24) 

and 

V e = V+ ^- a > '' A B F i) 2 ( 2 - 25 ) 



3 

i=i 



(l-a)5r 2 X " 
6 

i=i 
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respectively, with V the potential operator, a an arbitrary constant in the range 
[0, 1], A = h 2 /2m and Fj = VjV '. The resulting imaginary time propagator is accu- 
rate to order <5r 4 , and has been successfully applied to liquid 4 He in two and three 
dimensions. [T4"] This approximation offers also the advantage that adjusting the pa- 
rameter a it is possible to optimize the convergence, and a standard choice for a 
quantum system is a = 0.[H] A strategy to obtain a simpler, but equally accurate, 
approximation consists in applying a pair product assumption. [16J For sufficiently 
short time steps, in fact, the many-body propagator (in imaginary time) is well ap- 
proximated by the product of two-body propagators. [2] In this approximation, the 
small time propagator reads 



G(R m , Rm+i, St) =(4vrA(5r) 



-3N/2 



N 



Il eXP 



8=1 



ex P (~ u ( r ij,m, r ij,m+l)) 



X 



(2.26) 



where u is given as 



"('"mi 7"m+l) 



St 



v e (r m ) + 2v c (r m+ x)] m odd 
[2v c (r m ) + v e (r m+l )} m even. 



(2.27) 



The potentials v e (r) and v c {r) are defined as 



2 Q / dV 
v e {r) = V(r) + a-5r 2 X f — 

v c (r) = V(r) + (l-a) 1 -5r 2 X {^j 



(2.28) 



where V(r) is the potential experienced by two 4 He atoms at a distance r. The 
advantage is that there is no need to calculate Fi. As for the full Suzuki-Chin 
approximation, [13] also for the Pair-Suzuki the operators corresponding to physical 
observables must be inserted only on odd time slices in the imaginary time path. 

In order to fix the optimal small imaginary time step value, we have performed 
PIGS simulations with different initial wave functions. By considering decreasing 5t 
values with a fixed total projection time, r, we have taken the energy per particle 
E{t) as observable of reference. As an example, our results for SWF and CWF in 
the liquid phase are plotted in Fig. 2.2 We choose as optimal value 5r = 1/640 Kr 1 ; 



in fact, further reductions do not change the energy in a detectable way, i.e. within 



the statistical uncertainty. In Fig. 2J2 SWF and CWF do not converge to the same 
value simply because the considered total projection time r in this test is not enough 



to ensure convergence of E{r) to the ground state energy for CWF (see Fig. 2.3). 
Similarly, in the solid phase we take 5t = 1/960 Kr 1 . 
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Figure 2.2: Energy per 4 He atom E(t) vs. imaginary time step St. The total 
projection time is r = 0.1 K _1 . The calculations were carried out by projecting a SWF 
and a CWF for a system of 64 particles at the equilibrium density p = 0.0218 A~ 3 . 
Dashed lines are quartic fits to the data. Error bars, when not shown, are smaller 
than the used symbols. 



Once set the optimal 5t value, we have computed the diagonal properties of the 
system for increasing total projection time r until we reached convergence to a value 
that corresponds to the exact ground state result both for the liquid and for the solid 
phase. In the liquid phase we have computed also the one-body density matrix. 



PIGS results without importance sampling For the liquid phase we have pro- 
jected a SWF and a CWF. The energy per particle as a function of the total projection 



time r for both the wave functions is plotted in Fig. 2.3 We find that the energy 



converges, independently from the considered initial wave function, to the same value 
E = —7.17 ±0.02 K. This value, in spite of the small size of the considered system, is 
close to the experimental [T7] result E = —7.14 K. SWF converges very quickly, in fact 
r = 0.05 Kr 1 is already enough to ensure convergence. CWF instead, requires a three 
times larger imaginary time, i.e. r = 0.15 K _1 . Nevertheless, the quick convergence 
of also CWF is a really remarkable result. In fact, this means that PIGS efficiently 
includes the exact interparticle correlations through the imaginary time projections, 
without any need of importance sampling. Then, the choice of a good wave function, 
within the PIGS method, becomes a matter of convenience rather than of principle, 
since better initial wave functions only allow for a smaller total projection time r, 
and thus less CPU consuming simulations. 
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Figure 2.3: Energy per particle E as a function of the total projection time r obtained 
from PIGS simulations for liquid 4 He at the equilibrium density p = 0.0218 A~ 3 
by projecting a SWF (filled circles) and a CWF (open circles) and a GWF (open 
diamonds), r = result (filled circle) corresponds to the SWF variational estimate 
of E, the r = for the GWF is E = 122.08 ± 0.06 K and for CWF E is essentially 
infinite. Error bars are smaller than the used symbols. Dotted line indicates the 
convergence value E = — 7.17 ± 0.02 K. 



This convergence is confirmed also by the radial distribution function g(r) and 
the static structure factor S(k). For such quantities, the convergence rate is found to 
be similar to the energy one. In Fig. 2A we report the radial distribution function 
g(r) obtained by projecting both a SWF and a CWF at different imaginary time 
values. For r > 0.05 K _1 , SWF results at different r are indistinguishable within 
the statistical uncertainty (see Fig. 2.4 1). In fact, with SWF the exact result is 
reached within very few projection steps and then it is no more affected by further 
projections. As already pointed out, also CWF displays a fast convergence, as shown 
in Fig. 2A2, where Ag T (r) = gswF( r ) — PcwfC 7 ") * s shown. For increasing r, Ag T 
evolves toward a flat function, meaning that the systems described starting from the 
two different wave functions, i.e the strongly correlated quantum liquid of SWF and 
the ideal gas of CWF, are evolving into the same quantum liquid, which is the best 
reachable representation of the exact ground state of the simulated system. The same 
conclusion is inferred from the evolution of the static structure factor S(k), which is 



plotted in Fig. 2.5 
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Figure 2.4: Radial distribution function g{r) for bulk liquid 4 He computed in a 
cubic box with iV = 64 at the density p = 0.0218 A~ 3 with the PIGS method, a) g(r) 
obtained by projecting a SWF for r = 0.00, 0.05 and 0.25 K _1 . The r = 0.00 result 
corresponds to the variational SWF estimate of g(r). b) g(r) obtained by projecting 
a SWF for r = 0.25 K" 1 and a CWF for r = 0.25 K -1 . In the inset a zoom of the 
first maximum region, c) Ag T (r) = (?g WF (r) — S^wf! 7 ") at different r values, where 
(7g WF (r) is the g(r) computed by projecting a SWF for an imaginary time equal to 
r, and gcwF( r ) i s ^ ne same but by projecting a CWF. Note the smaller scale on the 
vertical axis 
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Figure 2.5: Static structure factor S(k) for bulk liquid 4 He computed in a cubic box 
with iV = 64 at the density p = 0.0218 A" 3 with the PIGS method, a) S(k) obtained 
by projecting a SWF and a CWF for r = 0.05 K -1 . b) S(k) obtained by projecting 
a SWF and a CWF for r = 0.40 K" 1 . c) AS T (k) = S^ WF (k) - S£ WF (k) at different r 
values, where Sg WF (k) is the S(k) computed by projecting a SWF for an imaginary 
time equal to r, and SQ WF (k) is the same but by projecting a CWF. Note the smaller 
scale on the vertical axis. 



PIGS results from a "bad" initial function In order to put a more stringent 
check on the PIGS method ability to converge to the exact ground state without any 
variational bias, we have considered also a "bad" initial wave function by projecting a 
GWF. Thus at the starting point of the imaginary time path there is now a strongly 
localized Einstein crystal. We note that, differently from the other considered cases, 
the GWF is not Bose symmetric; as consequence of this choice, the projection relation 
(2.8) is not Bose symmetric and thus requires symmetrization. The symmetrization 
has been introduced with the sampling of the permutations between polymers; this 
is a standard technique used in Path Integral simulations and will be described in 
Sec. EH 

With this initial function, we find even in this case that the energy converges to 
the same value as before (see Fig. 2.3). Thus PIGS is able not only to drop from the 
initial wave function the wrong information of localization, but also to generate at 
the same time the correct correlations among the particles. GWF needs r = 0.5 K _1 
to converge, which is ten times larger than the SWF value. 
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Figure 2.6: Radial distribution function g{r) for bulk liquid 4 He computed in a 
cubic box with iV = 64 at the density p = 0.0218 A -3 with the PIGS method, a) 
g(r) obtained by projecting a SWF for r = 0.40 K -1 and a GWF for r = 0.50 K" 1 . 
In the inset a zoom of the first maximum region, b) Ag T (r) = gswp{r) — PcwfC 7 *) a ^ 
different r values, where ^swfM i s the g(r) computed by projecting a SWF for an 
imaginary time equal to r, and gQ WF (r) is the same but by projecting a GWF. Note 
the smaller scale on the vertical axis. 



Again this convergence is confirmed also by the radial distribution function g(r) 
and the static structure factor S(k). In Fig. 2.6 we report the radial distribution 
function g(r) obtained by projecting a GWF at different imaginary time values com- 
pared with the ones coming from the projection of SWF. It is evident that small 
imaginary time is not enough to leave out the wrong information in the GWF. For 
lower t values, there are still reminiscences of the starting harmonic solid, which are 



progressively lost as the projection time increases. This is made clearer in Fig. |2.6b 
where we plot the difference Ag T (r), at fixed imaginary time r, between the g(r) 
computed by projecting the SWF and the one obtained by projecting the GWF. A 
similar behavior is observed in the evolution static structure factor S(k), plotted in 
Fig. 2.7 For the GWF, the Bragg peak shown at small r values (Fig. |2. T a,) , which is 
typical of the solid phase, becomes lower and lower as the projection time is increased 
(Fig. [2?7», until convergence is reached (see Fig. [277c 



From the plot of the energy per particle vs. the total imaginary time r it is 
possible to estimate the overlap per particle of the initial wave function on the exact 
ground state. [18] By using the results in Fig. 2.3 we find that the overlap of SWF 
is about 99%, while the GWF one is about 10%. That SWF has an high overlap 
with the ground state is not a surprise; it was qualitatively expected since SWF is 
presently the best available wave function for He. [6] However a 99% overlap is really 
remarkable and provides a further argument on the goodness of SWF. On the other 
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Figure 2.7: Static structure factor S(k) for bulk liquid 4 He computed in a cubic box 
with iV = 64 at the density p = 0.0218 A" 3 with the PIGS method, a) S(k) obtained 
by projecting a SWF and a GWF for r = 0.05 K _1 . It is evident in the GWF result 
the presence of the Bragg peak. Note the logarithmic scale, b) S(k) obtained by 
projecting a SWF for r = 0.40 K" 1 and a GWF for r = 0.50 K _1 . The Bragg peak is 
no more present in the GWF result, c) AS T (k) = Sg WF (k) — SQ WF (k) at different r 
values, where Sg WF (A;) is the S(k) computed by projecting a SWF for an imaginary 
time equal to r, and <Sq WP (A;) is the same but by projecting a GWF. Note the change 
of the vertical scale. Error bars are smaller than the used symbols. 
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Figure 2.8: One-body density matrix pi obtained from PIGS simulations for liquid 
4 He at the equilibrium density p = 0.0218 A -3 by projecting a SWF, a JWF and a 
GWF for an imaginary time r = 0.30, 0.40 and 0.80 K _1 respectively. The dotted 
line indicates the condensate value uq = 0.069 obtained from an independent PIGS 
simulation. [19] 

hand, a poor overlap of GWF was somehow expected, since the parameter C was 
chosen to strongly localize the atoms of the bulk liquid around fictitious equilibrium 
positions on a regular lattice. 

Off-diagonal properties Besides the diagonal ones, also off-diagonal properties, 
such as the one-body density matrix, are accessible within PIGS simulations. The 
one-body density matrix pi(f,r f ) represents the probability amplitude of destroying 
a particle in f and creating one in ? . Its Fourier transformation represents the 
momentum distribution. In first quantization pi is given by the overlap between 
the normalized many-body ground state wave functions ipo(R) an d ipo(R'), where the 
configuration R' = (f*, f 2 , ■ ■ ■ , rjv) differs from R = (r, r 2 , . . . , rjv) only by the position 
of one of the N atoms in the system. If if>o(R) is translationally invariant, p\ only 
depends on the difference \r — r*|, thus 

Pl ^~7^) = N J dr 2 ...dr N %(R)MR')- (2-29) 

The Bose- Einstein condensate fraction n is equal to the large distance limit of p%(f — 
f > ). In fact, if pi has a nonzero plateau at large distance, the so called off-diagonal 
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long-range order (ODLRO), its FT contains a Dirac delta function, which indicates a 
macroscopic occupation of a single momentum state, i.e. Bose-Einstein condensation. 



The exact p\ can be obtained in PIGS simulation by substituting ip in (2.29) 
with ip T with r large enough. This corresponds to the simulation of a system of N — 1 
linear polymers plus a polymer which is cut into two halves, called half-polymers, 
one departing from r and the other from r*. Thus pi is obtained by collecting the 
relative distances among the cut ends of the two half-polymers during the Monte 
Carlo sampling. The present computation of p\ has been obtained by implementing 
a zero temperature version of the worm algorithm. [2U] We have worked with a fixed 
number of particles and not in the grand canonical ensemble, similarly to what has 
been done at finite temperature in Ref. [IS]. In practice this corresponds to a usual 
PIGS calculation of p\ where "open" and "close" moves have been implemented [20] 
in order to visit diagonal and off-diagonal sectors within the same simulation. The 
advantage of doing this does not come from the efficiency of the worm algorithm to 
explore off-diagonal configurations, because similar efficiency is obtained with PIGS 
when "swap" moves are implemented. [11] The benefit in using a worm-like algorithm 
here instead comes from the automatic normalization of p\ which is a peculiarity of 



this method. |20j In Fig. 2.8 we report p\ obtained in PIGS simulations of bulk liquid 



He at p = 0.0218 A~ 3 by projecting either a SWF, a JWF and a GWF. All the 



simulations give the same result, shown in Fig. 2.8 which turns out to be compatible 



with the recent estimate obtained with PIGS given in Ref. [19] of no = 0.069 ± 0.005. 

Results on the solid system We have performed the computation at density 
p = 0.0313 A -3 , where 4 He is in the solid phase, by projecting a SWF and a CWF. 



Our results for the energy per particle are plotted in Fig. 2J3 as a function of r. In 
both cases we find convergence to the value E = —5.34 ± 0.02 K. Even in this phase 
the convergence of SWF is faster, being r = 0.05 K -1 enough to reach convergence. 
In the case of CWF convergence is reached only for a much larger imaginary time 
t = 0.80 K _1 . 

Also in this case convergence is obtained for the radial distribution function and 



for the static structure factor, reported in Fig. 2.10 and Fig. 2.11 respectively. From 



Fig. 2.10 1 it is evident that SWF has reached the true ground state with few projection 
steps, since the results for g(r) at r = 0.05 K _1 and r = 0.80 K are indistinguish- 
able. The evolution toward the correct ground state of the projected CWF is instead 
detectable. The presence of the crystalline structure is mainly evident in the static 
structure factor, where a Bragg peak grows with increasing r (see Fig.|2~Tl>,b). The 



emerging of the correct solid structure by projecting a really poor wave function such 
as the CWF is made evident by the trend toward a flat function of the differences 



Ag T (r) and AS T (k) plotted in Fig. 2.10: and Fig. 2.11: respectively. 



Conclusions of the test phase In this section we have studied with the Path 
Integral Ground State method diagonal and off-diagonal properties of a strongly in- 
teracting quantum Bose system like the bulk liquid and solid phases of 4 He. We have 
obtained convergence to the ground state values of quantities like the total energy, 
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Figure 2.9: Energy per particle £ as a function of the total projection time r 
obtained from PIGS simulations of an fee 4 He crystal at the density p = 0.0313 A~ 3 
by projecting a SWF (filled circles) and a CWF (open circles). Dashed line indicates 
the convergence value E = —5. 34 ±0.02 K. 

the radial distribution function, the static structure factor and the one-body density 
matrix projecting radically different wave functions: equivalent expectation values in 
the liquid phase have been obtained using as initial wave function a shadow wave 
function, a Gaussian wave function with strongly localized particles of an Einstein 
solid without interparticle correlations and also a constant wave function where all 
configurations of the particles are equally probable. Similarly in the solid phase equiv- 
alent expectation values have been obtained by considering a shadow wave function, 
which describes a solid, and a constant wave function which describes an ideal Bose 
gas. The present analysis demonstrates the absence of any variational bias in PIGS; 
a method that can be thus considered as unbiased as the finite temperature PIMC. 
This remarkable property comes from the accurate imaginary time propagators, ex- 
actly the same used with PIMC, that do not depend on the initial trial state. It 
remains true that the use of a good variational initial wave function greatly improves 
the rate of convergence to the exact results. 

We have addressed here only the case of a realistic interaction potential among 
Helium atoms. However one can reasonably expect that this conclusion holds even for 
very different kinds of interaction, once an accurate approximation for the imaginary 
time propagator is known (for example hard- spheres (21] or hydrogen plasma[22j). As 
far as pathological potentials like the attractive Coulomb one are concerned, PIGS 
would suffer the same limitations of PIMC if inaccurate approximations of the prop- 
agator were used. [23] 
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r (A) 

Figure 2.10: Radial distribution function g{r) for bulk solid 4 He computed in a 
cubic box with N = 32 at the density p = 0.0313 A -3 with the PIGS method, a) 
g(r) obtained by projecting a SWF for r = 0.05 and 0.80 K -1 . b) g{r) obtained 
by projecting a SWF and a CWF for r = 0.80 K _1 . In the inset a zoom of the 
first maximum region, c) Ag T (r) = <?g WF (V) — ^ WF (r) at different r values, where 
gg WF (r) is the g(r) computed by projecting a SWF for an imaginary time equal to r, 
and gcwF( r ) * s the same but by projecting a CWF. 

2.2 Path Integral at finite temperature 

Up to now we have focused on the problem of evaluating T = K expectation 
values; the formalism of the previous section, however, can be used with very small 
modifications also for quantum thermal averages, the resulting methodology is named 
Path Integral Monte Carlo (PIMC). This methodology was developed well before 
PIGS in the work in Ref. |24J. The physical properties of the system are obtained 
from the thermal density matrix 

p = — (2.30) 
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Figure 2.11: Static structure factor S{k) for bulk solid 4 He computed in a cubic 
box with iV = 32 at the density p = 0.0313 A" 3 with the PIGS method, a) S(k) 
obtained by projecting a SWF and a CWF for r = 0.05 K _1 . b) S(k) obtained by 
projecting a SWF and a CWF for r = 0.80 K . The black dots are under the red 
ones, c) AS T (k) = S^ WF (k) — SQ WF (k) at different r values, where Sg WF (k) is the 
S(k) computed by projecting a SWF for an imaginary time equal to r, and SQ WF (k) 
is the same but by projecting a CWF. Error bars are smaller than the used symbols. 
Notice the logarithmic scale in panels a) and b). 
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where (3 = fc& is the Boltzmann constant and the normalization Z = Tr (p) is 
the partition function of the system. The expectation value of an observable O is 



< O >-- 



Tr Op 



(2.31) 



It is evident in Eq. (2.30) that the unnormalized density matrix operator is formally 



identical to the quantum imaginary time operator appearing in Eq. (2.5) if one 
chooses t = p. The density matrix p in coordinate representation becomes 



G(R,R',p) = (R 



R! 



(2.32) 



Fixed this set of basis, |-R)(-R| = / dR, the trace operator acts on the density matrix 
as follows 



Tr (p) = J dRG {R, R,f3) = Z 
1 



(0) = — I dRO (R) G (R, R, (3) 



(2.33) 
(2.34) 



where O is diagonal in the coordinate representation. 

Using the Path Integral notation and following the same procedure employed for 



PIGS, Eq. (2.31) becomes 



O 



1 

Z 



/M M-l 
11 dR t O (R k ) J] G (Rj,R j+ i, St) G (R m , R x , St) 
i=i j=i 



(2.35) 



where here St = (3/M; the cyclic property of the trace operation allows O(R) to be 
evaluated at any position in the path integral, in other words 1 < k < M. 



Eq. (2.18), however, does not posses the Bose symmetry; in order to introduce 



either the Bose or the Fermi statistics one has to symmetrize the density matrix (2.35) 
over the permutations of the particle labels. 



G B (R, R',P) = ^ J2 g pR > P) 

p 

G F (R, R', 0) = Yl G ( R '> PR > P 



(2.36) 
(2.37) 



the first equation holds for the Bose statistics and the second for the Fermi statistics. 
The permutation operator P acts on the many-body coordinate R = {fi]^ =l by 

applying a cycle of n P exchanges between particle indices, PR = {rp(i)} =1 - The 
sum in the two equations is meant as a sum over all the iV! possible permutations. 

The Fermi symmetrization introduces negative density matrices that no longer 
can be directly interpreted as probability densities; there are different techniques [25] 
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that can re-express the density matrix as a definite positive object with a weight on 
the sampled configurations that introduces the sign given by the Fermi symmetry, 
however these techniques yield a signal to noise ratio [2Sj that dramatically decreases 
as e~ lN ; as consequence of this sign problem, these methodologies are restricted to 
small particle numbers. 



Due to its affinity with Eq. (2.11 ), Z is also the partition function of a system of 
classical closed polymers. 




Figure 2.12: Representation of the polymers in PIMC and the effect of a permutation; 
a permutation between two polymers results in a new compound polymer with twice 
the length in imaginary time. This is a configuration that can not be sampled with 
moves that involve only a single polymer. Grey beads and lines represent the removed 
segment of the polymer. This picture can also be viewed from right to left: in this 
case, grey beads and lines represent the new segment of the polymer; this permutation 
splits a ring polymer of length 2/5 into two ring polymers of length /3. 



The differences between PIGS and PIMC in "polymer language" are minimal; first, 
in PIMC does not appear any trial wave function and the quantum imaginary-time 
is proportional to (3 = l/(fc^T); a more subtle difference that has deep consequences, 
however, is the different topology of the polymers. Due to the trace operation, in 
fact, the polymers in PIMC are represented by Eq. (2.33); these polymers close on 
themselves (for this reason they are also called ring polymers, see Fig. 2.12) and the 
role of permutations becomes more important than in the PIGS case. This is so be- 
cause Eq. (2.33) is not yet Bose-symmetrized and the effect of the symmetrization 
(2.36) is that the polymers representing the particles of the system no longer close on 
themselves but, instead, are allowed to close on another ring polymer; this is explicitly 
shown in Fig. 2.12 The relevance of the symmetrization in PIMC is also noted by 



the fact that permutations explore topologically different configurations of the system 
that could not be obtained with single polymer sampling. This situation is different 
from the PIGS case; in PIGS, in fact, permutations between polymers will simply 
yield another configuration of open polymers that is topologically identical and can 
be obtained with single polymer moves; this implies that the sampling of permu- 
tations is not necessary if the Bose symmetry is already introduced by the chosen 
trial wave function. Another difference is that open polymers have less constraints 
on their structure and, compared to ring polymers, can be moved more efficiently by 
the Metropolis sampling; this is especially true when the probability distribution to 
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sample has many local maxima: in this case, with open polymers is generally easier 
to satisfy the ergodicity of the sampling. 
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Chapter 3 



Polarization energy of 
two— dimensional 3 He 



In this chapter, the subject under study is a system of two-dimensional 3 He (2d 3 He) 
in a wide range of densities in the liquid region and up to freezing. This subject 
is of interest because, as it has been shown in Ref. p], a pure 2d 3 He system is a 
good approximation of a quasi-two-dimensional 3 He sample. Such system can be 
experimentally realized over a wide range of liquid densities by adsorbing 3 He on a 
variety of preplated graphite substrates [21 El H]. Regarding the effective mass m* 
and the spin susceptibility x/Xoi the system behaves in good approximation like a 
perfect Fermi liquid pQ: the enhancement of x/Xo increases with the density and m* 
is consistent with a divergence near freezing density. This behavior at freezing has 
been interpreted [1] as a signature of a Mott transition leading to an insulating crystal. 
Theoretical studies [5], however, suggest that the singularity of m* and freezing could 
not have the same origin and the freezing density is influenced by the preplated 
substrate. In this context, it is relevant the study of the effect of correlations without 
any effect induced by the external potential of the substrate. 

Bulk 2d 3 He is interesting also from the theoretical point of view because, being a 
strongly interacting system at high densities, it provides a severe test case for micro- 
scopic calculations [5]. Some of the most powerful tools to study strongly interacting 
systems are QMC methods. The so-called fixed-node (FN) approximation [7j, used 
in most Fermionic QMC calculations, however, has been argued to give a significant 
bias in the polarization energy of three-dimensional liquid 3 He[8] at high density. 
We have thus performed QMC simulations beyond the FN level, following a formally 
exact method^] that is referred here as Fermionic Correlations (FC). With both the 
FN and FC methods we have calculated the ground-state energy per particle e = j| 
of the 2d 3 He liquid at zero temperature as a function of the number density p and 
the spin polarization £. 

The FC method is slightly different from the well known transient estimate (TE) 
technique |1U], the basic idea is to perform simulations relying on the basic Hamil- 
tonian in an enlarged, unphysical space of states of any symmetry, including those 
with Fermi and Bose statistics. The ground state energy of the physical fermionic 
3 He is considered as an excitation energy of the absolute bosonic ground state, which 
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is sampled exactly with QMC. In this approach one 'trades' (in a sense which we will 
explain below) the sign problem faced by TE[T0] for the analytic continuation needed 
to extract excitation energies from suitable imaginary-time correlation functions. A 
mixed approach, devised to ease detection of the asymptotic convergence of TE by a 
Bayesian analysis of imaginary-time correlation functions, was proposed By Caffarel 
and Ceperlev[TT]. 

We compared our results also with a previous FN QMC calculation [12] that how- 
ever is limited to low densities and only considers the paramagnetic fluid phase. We 
find indeed that the FN level of the theory and the exact calculation predict a qual- 
itatively different behavior. This is rather expected because the accuracy of the FN 
approximation in the high density regime is questionable [8]. In fact, within FN the 
system becomes ferromagnetic well before crystallization takes place upon increasing 
the density, whereas the unbiased calculation shows that the spin polarization of the 
fluid is preempted by freezing, as observed experimentally. From the estimated curve 
e(Q we obtain a spin susceptibility enhancement in quantitative agreement with the 
available measurements. 

3.0.1 QMC simulation 

We simulate iV particles with the mass 7713 of 3 He atoms, interacting with the HFDHE2 
pair potential[13] in periodic boundary conditions, which is the most accurate two- 
body potential for Helium systems[14j. The simulation box, of area Q, is a square 
of side L for the liquid phase; for the solid it is a rectangle which accommodates a 
triangular lattice. The Hamiltonian is 

fc2 N N 
3 i=l i<j=l 

The simulations of the bosonic 3 He were made with the SPIGS technique, using 
the Pair Product approximation for the propagator at short imaginary times. For 
this system, we have verified that, using a Shadow Wave Function, a projection time 
of t = 0.2 K^ 1 is enough to yield an accurate description of the ground state. As for 
the time-step, we used 5t = 1/160 K -1 and verified that it is small enough to be an 
accurate approximation. 

Fixed— Node approach 

The fixed-node approximation [7J is one of the most commonly used approaches in 
the QMC simulation of Fermi systems. FN stochastically solves the imaginary-time 
Schrodinger equation subject to the boundary conditions implied by the nodal struc- 
ture of a given trial function ty T . This approach gives a rigorous upper bound to the 
ground state energy, which often turns out to be extremely accurate. 
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For a given spin polarization, i.e. considering iV-j- spin-up and N± = N — Nf 
spin-down 3 He atoms, is chosen of the form 



* T (1l) = V{T^)V{n^j{n) X< i (3-2) 

where 1Z = (fi, ■■■,r N ), TZ^ = (ri, ...,fjv t ), ^4. = (f/v t +i, •••,^v), and the whole depen- 
dence on the spin degrees of freedom is contained in \Qi a s P m eigenfunction for the 
given polarization 

C = ^^ , (3.3) 

The Jastrow factor, 

= II ex P (-^(l^-^l)) . (3-4) 

i<j \ / 

describes pair correlations arising from the interaction potential; we use a simple 
pseudopotential of the McMillan form u(r) = (b/r) 5 . Finally, the simplest form of 
the antisymmetric factors T> (7fyg,) is in the form of Slater Determinants of plane 
waves: 

V (TZ U ) = det N exp(iki ■ ^)}. J (3-5) 

More accuracy in the FN results is achieved by introducing also backflow correlations 
via quasi-particles vector positions: 



V (Tl u ) = det ^{exp(i£ • ^)}. J (3-6) 

^ = ?j + £^ i=1 V(\rj ~ n\) (fjf - . 
For the backflow correlation function 77 (r) we adopt the simple form: 

77(7-) = Ae~ Bir - c)2 . (3.7) 

We will refer to the two choices respectively as plane waves fixed-node (PW-FN) 
and backflow fixed-node (BF-FN). For each density, the variational parameters b, A, 
B and C are optimized using correlated sampling [16] at ( = 0, and left unchanged 
at different polarizations. The backflow parameters, for each density, are shown in 
Table O 

Part of the bias related to the finite size of the simulated system arises from shell 
effects in the filling of single-particle orbitalsfTT]. This bias can be substantially 



reduced adopting twisted boundary conditions [T7]. i.e. choosing k appearing in (3.5) 



and (3.6) inside the set: 

-> 2vm + 9 



l (3 ' f 
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Table 3.1: Backflow parameters used for each density 



Density ( A -2 


A 


B ( A' 1 ) 


C(A) 


0.020 


1.90393 


0.117865 


-1.89877 


0.050 


1.124523 


0.112559 


-0.94888 


0.060 


1.017654 


0.147372 


-0.51614 



where n is an integer vector, L is the side of the simulation box f2 and 6 is a twist 
parameter Q{ G [0, it] which, at the end of the calculations, is averaged over. 

In the solid phase, quantum exchanges are strongly suppressed and the energy 
difference between a Fermionic and a Bosonic crystal is negligibly small for the pur- 
pose of locating the liquid-solid transition. We will therefore replace the energy of 
3 He with that of a fictitious bosonic Helium of mass m.3, which can be calculated 
exactly [TBT EH 120]. The small error incurred by such replacement is bound by the 
difference between the fermionic Fixed-Node (FN) energy and the unbiased bosonic 
energy. This difference, calculated[21J as a check at the melting density where it is 
expected to be largest, is indeed in the sub-milliKelvin range. 

We stress that we actually made a particular choice of trial wave functions; the 
obtained results depend on such a choice: when we will speak about 'fixed-node 
level' of the theory or about 'fixed-node approximation' we will always implicitly 
refer to the above mentioned trial wave functions. Naturally it could be possible to 
improve the fixed-node results using more sophisticated wave functions; instead, we 
have chosen to follow another way with the FC method; this method is in principle 
exact and does not depend on a particular choice of the wave function. 



Fermionic correlations approach 

As mentioned before, for the fluid phases the FN approximation may not be accurate 
enough, particularly at high density where correlations are stronger and the energy 
balance between different polarization states is more delicate. 

In order to go beyond the FN level and obtain accurate data, we use the FC 
technique [9] which is in principle exact, even if limited to small system sizes. 

The idea, with similarities with the transient estimate formalism [TOT ITT] , is that 
of viewing (3.1) as an operator acting inside the Hilbert space f-L(N) = (L 2 (f2)) 07V , 
that has no constrains on spin and statistics: one can use Quantum Monte Carlo to 
sample the lowest energy eigenfunction ipo(TZ) of H among the states of both Bose 
and Fermi symmetry. 

It is known [22J that ipo must share the Bose symmetry of the Hamiltonian, so 
that: 

B (%\Hj) ) n{N) 

is the Ground State energy of a fictitious system of N Bosons of mass interacting 
via the potential v(r). 
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The connection between the fermionic energies is retrieved in the following way: 
let us fix a spin polarization, it is surely a good quantum number since the basic 
Hamiltonian is spin-independent. As discussed in Ref. jS], if we are able to define an 
operator Af such that, inside T-L(N), 

iPf(K)=(AfiJ>o)(TI) (3.10) 

has non-zero overlap with the configurational part of any exact fermionic Ground 
State of H for the given (, we can define the imaginary-time correlation function: 

(Vol (e T6 AU- TA )AFi>o)H(N) 
which can be straightforwardly evaluated in a bosonic QMC simulation 



This is readly made because the evaluation of Eq. (3.11) at a certain discrete imagi- 
nary time 77 = ISr can be done with the evaluation of the Slater determinant Af on 
two time sectors located at different imaginary times tq = mbr and t% = (m + l)8r. 
The actual calculation has been done using the Path Integral Ground State with 
Shadow Wave Functions (SPIGS) that has been described in chapter [5J The lowest 
energy contribution in Cf(t) provides the exact gap between the fermionic and the 
bosonic ground states of the two-dimensional Fermi liquid; this can be readily seen 



by formally expressing (3.11 ) on the basis {il>n}n>o of eigenvectors of H corresponding 
to the eigenvalues {E n } n > : 



- r (^-^B\\(A F 'ipo\^n}n(N)\ 



n=0 



Cf(t) = 2^ e- r ^"-^M , ( 3 - 12 ) 



A quite natural choice [S] is to define Af borrowing suggestions from the form of the 
trial wave function for the FN calculation, i.e.: 

(Xvo) (n) = f v(n t )v(n x )M^) (3.13) 



where we can choose either the definition (3.5) of T> or the definition (3.6). We will 



refer to such choices simply as the plane waves fermionic correlations (PW-FC) and 
the backflow fermionic correlations (BF-FC). Naturally the final results for the Bose- 
Fermi gap should coincide within statistical uncertainties, and the actual comparison 
can be a good test for the robustness of the approach. 

We observe that the sign problem is not really avoided as it manifests itself in two 
ways: on one hand poor choices of the wave functions appearing in the correlation 
functions imply the necessity to consider very large r regions of the correlation func- 
tion; on the other hand, since the gap energy is an extensive quantity, the exponential 
decay of the correlation functions increases in the limit N — > 00, making impractical 
the extraction of meaningful information. 
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3.0.2 The Bosonic System 



Figure |3.1| shows the state equation of both the solid and the liquid phases of the 
system; in Table 3.2 the values of the energies at each density are shown. 



Table 3.2: Potential (E pat ), Kinetic (E kin ) and Total (E tot ) Energy per particle of 
"bosonic" 2d 3 He at each studied density. The system has N = 26 atoms in a square 
box of late L = (N/p) 1 ^ 2 . Tail corrections to the potential energy to account for the 
finite size of the system have been applied only to E to t/N; instead, for a more direct 
comparison, E pot /N is exactly the output of the simulations. 



Density ( A~ 2 


E pot /N ( K ) 


E kin /N ( K ) 


E tot /N ( K ) 


0.015 


-1.176 (5) 


1.107 (7) 


-0.10 (1) 


0.020 


-1.640 (6) 


1.58 (1) 


-0.10 (2) 


0.025 


-2.150 (6) 


2.12 (1) 


-0.09 (2) 


0.030 


2.75 (2) 


-2.709 (9) 


-0.04 (2) 


0.035 


3.47 (2) 


-3.318 (7) 


0.06 (4) 


0.040 


4.29 (2) 


-3.97 (1) 


0.22 (4) 


0.045 


5.26 (2) 


-4.71 (1) 


0.45 (4) 


0.050 


6.35 (2) 


-5.47 (1) 


0.78 (4) 


0.055 


7.62 (3) 


-6.28 (1) 


1.25 (4) 


0.060 


9.16 (3) 


-7.12 (1) 


1.86 (4) 


0.065 


10.81 (3) 


-7.98 (2) 


2.67 (5) 


0.070 


12.75 (4) 


-8.91 (2) 


3.66 (7) 


0.075 


15.00 (3) 


-9.94 (3) 


4.77 (6) 


0.080 


17.27 (4) 


-10.84 (2) 


6.12 (7) 


0.085 


19.70 (4) 


-11.65 (2) 


7.81 (7) 


0.090 


22.32 (5) 


-12.28 (3) 


9.90 (9) 


0.095 


25.19 (6) 


-12.72 (3) 


12.5 (1) 


0.100 


28.31 (7) 


-12.88 (3) 


15.7 (1) 



The agreement with data in Ref. [12] is good (Fig. 3.2), although the comparison 
must take into account the slight difference coming from the interaction potential used 
in our simulation (Aziz '79, Ref. p2]), and that used in Ref. p2] (Aziz '87, Ref. [21]). 
This difference, E$ 7 — E 79 , depends on the density p and can be estimated with 
the following computation if one neglects the dependence of the radial distribution 
function, g (r), on the interaction potential v (r): 



CO 



E S7 - E 79 -np dr [v 87 (r) - v 79 (r)] g (r) (3.14) 
J o 

With the Maxwell construction obtained from the polynomial fit of the equation 



of state in the liquid and the solid phase (see Table 3.3 for the fitting parameters), 
the freezing point is estimated at a density of 0.069A~ 2 , while the melting point is 
approximately at a density of 0.073A -2 . 
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Figure 3.1: State equation of the bosonic system used to calculate the correlation 
functions. The lines represent a fit of the data in both the solid (red line) and the 
liquid (black line) phase. The error bars are hidden inside the symbols. These results 
are in good agreement with Ref. [T2] 



Table 3.3: Fitting parameters for the solid and the liquid phases. The polynomial 
that has been fitted to the data is of the form E/N = ap 5 + bp A + cp 3 + dp 2 + ep. 



KA^ 



b ( KA i 



KA^ 



d ( KA 4 



Phase 



a 



c 



d ( KA 2 



liquid 
solid 



1745632.16048 
3063129.8758 



-66748.2524 
-532201.7865 



9553.61311 
58695.29501 



-61.70852 
-2598.78493 



-7.51343 
57.3547 



High accuracy in these computations is very important for two reasons: first, 
even though the bosonic system is unphysical in itself, its energy is required for the 
evaluation of the energy of the real 3 He system, and finally, high quality of data is 
essential for a successful inversion of the Laplace Transform. 

3.0.3 The Twist Averages 

One of the main differences between the TE and the FC methods is the way in which 
twisted boundary conditions (TBC) are used. In FC, TBC are not applied whenever 
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Figure 3.2: Comparison between the equation of state of mass-3 bosons obtained 
with the Diffusion Monte Carlo and that obtained with the SPIGS. The systematic 
difference of roughly 0.03K is due to the different interaction potential that has been 
employed. 



a particle crosses a boundary of the simulation box but are taken into account during 
the evaluation of the correlation function; more precisely, the twist angle is introduced 
during the preparation of the Slater matrices, namely: 



S 



Ui4 


L J 




if fci-4 


4) 








e V 





\ 



i \ A - : iv t + 7 L )-f N . 



(3.15) 



7 



where Q\ is a given twist angle and is the number of spin up particles; naturally 
the same applies for the spin down particles. The coordinates f, can be, like in the 
FN case, with or without backflow correlations, in the latter case we place fj = r z abs , 
where the superscript means that the coordinates of the i-th particle are obtained 
without invoking the periodic boundary conditions (pbc); this is done in the following 
way: each particle has two coordinate systems, the first are coordinates inside the 
simulation box, which we refer here with the superscript "pbc" , these are obtained by 
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invoking the pbc whenever the particle crosses a boundary of the simulation box, the 
second coordinate system is made with absolute coordinates that are not constrained 
in the simulation cell: if a particle moves from a position xf bc to a position x^ ° 
its absolute coordinate changes accordingly x 
last subscript means that the displacement (x 



abs 
pbc 



X 



abs 



+ [X 



pbc 



— X 



pbc\ 



%i c )pbci where the 
is calculated with periodic 



boundary conditions. The other choice for fj is with backflow correlations, 



A? 



r . 
1 i 



abs 



+ E * 

j=l;jy*=i 



Tn —Ti\, 

J L I pbc 



'pbc 



(3.16) 



where r\ (r) is defined in Eq. (3.7). 



Another delicate point in the construction of the Slater matrix (3.15) is the choice 



of the values {k n + for Q\ = the choice would reduce to the wave vectors 

{k n } inside the Fermi surface corresponding to iVf particles, however it may not be 
the case for certain choices of the twist angle 61, the procedure to follow is that, for a 
given twist angle 9i, the wave vectors are those that give the first iV-f lowest energies 



E n = A 



01 
L 



(3.17) 



where A = h 2 /(2m3 He ). 

For the evaluation of a single energy gap, 15 different correlation functions have 
been used for every Monte Carlo block. Each correlation function corresponds to a 
twist angle. This choice leads to a uniform distribution of of the twist angles in an 
area of the first Brillouin zone of the simulation box that contains no symmetries. 



This area is shown in Fig. 3.0.3 



Following the prescription in Eq. (3.13), once the Slater matrices have been pre- 



pared, the Slater determinant has to be computed. This operation will be invoked 
many times during a Monte Carlo step and an efficient algorithm is advised. Our 
choice has been the LU decomposition that writes a matrix A as a product of an 
upper triangular matrix U and a lower triangular matrix L. The determinant of A 
is then the product of the diagonal elements of U and L. The LU decomposition 
is extensively described in Ref. [25] for real matrices, the case of complex matrices 
is easily generalizable: from the algorithm in Ref. [25] it is enough to redefine the 
matrices L and U to complex matrices L = $IL + and U = $tU + QU; all the 
algebraic operation have then to be ambiented in the complex field and the resulting 
determinant will be a complex number. A complex number for the Slater determinant 



will yield a complex imaginary time correlation function; however, due to Eq. (3.12) 



this correlation function will have, on average, an imaginary part compatible with 
zero. 

The procedure to obtain the polarization curves may be schematized in the fol- 
lowing steps: 
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Figure 3.3: Schematic representation of the chosen twist angles in the first Brillouin 
zone of the simulation box. The simulation box is a square of late L. 



• Consider all of the correlation functions corresponding to a given twist angle 
that have been computed for each block of the simulation. Calculate the error 
using the central limit theorem. 

• From the previously calculated error, infer the error of the correlation function 
relative to a single block. 

• Apply the inversion method and localize the position of the first peak. In this 
work we used the GIFT method but we obtained compatible results also with a 
fit of the long imaginary-time part of the imaginary-time correlation function. 

• Repeat the previous steps for every twist angle and hence perform a weighted 
average according with the symmetries of the first Brillouin zone. This yields a 
single-block estimate of the energy gap. 

• The final energy-gap value for a given polarization is the block average of the 
single-block estimates that have thus been obtained. 

The described procedure has been applied for every polarization. This implies 
thousands of Laplace transforms to be inverted with relative peaks to be localized; it 
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Figure 3.4: (Color online) Upper panel: Imaginary time correlation functions, Cp(r) 



corresponding to the two different choices of the operator in (3.13). Lower panel: 



reconstructed spectral functions sp(u) obtained with the GIFT method. 



is a task that is unlikely to be done manually and an automated procedure was imple- 
mented. However, particular attention must be paid in the automatic localization of 
the peaks in order to avoid false values and thus systematic errors. Possibly, different 
localization algorithms must be applied and a check by eye of a random selection 
of the inversion results is highly advised. We would like to emphasize again that 
the most delicate part of this method lies in the numerical inversion of the Laplace 
transform of the correlation functions but also in the data analisys of the results. 

3.0.4 Analytic Continuation 

Once we have achieved a QMC evaluation of Cf(t), the information about the Bose- 
Fermi gap Abf = Eq — Eq is contained in the resulting correlation functions. The 
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results for Cf{t) appear as simple smooth decreasing functions, whose values can be 
evaluated only in correspondence with a finite number of imaginary-time values, say 
{t , Ti, r 2 , Ti}] moreover the data are perturbed by unavoidable statistical uncer- 
tainties. The Bose-Fermi gap Abf is thus hidden inside the sets of limited and noisy 
data. How can we extract it? 



In the upper panel of Fig. |3.4| we show two imaginary time correlation functions 
Cf(t), respectively a PW-FC and a BF-FC, corresponding to the same spin polar- 
ization and twist parameter. The long-r tails of the two curves tend towards a linear 
behavior (in logarithmic scale) with the same slope, and this is a general feature 
shared by all the functions we have evaluated. This indicates that, because of the 
finite-size of the system (and selection rules on the total momentum), the fermionic 
spectrum has a sufficiently defined gap, i.e. the lowest-energy term exp(— Abft) 



in the correlation function (3.12) appears to be quite well resolved with respect to 
contributions from higher fermionic energies. The difference between the two curves 
(in particular the rigid shift between their asymptotic tails) arises from the spectral 
weight of the Ground State contribution, which is higher when backflow correlations 
are taken into account, as expected. 

In this favorable situation, the Bose-Fermi gap can be reliably extracted by simply 
fitting an exponential to the long-time tail of the correlation function. 

This key result is strongly supported by a more sophisticated approach, which 
evaluates Abf by performing the full Laplace transform inversion of Cf{t), i.e. solving 

r+oo 

C f {t) = / due- TUJ s F {uj) , (3.18) 
Jo 

for the unknown spectral function sf{oo)- Recently a new method, the genetic inver- 
sion via falsification of theories (GIFT) method[26j, has been developed to face general 
inverse problems and in particular it has allowed to reconstruct the excitation spec- 
trum of superfluid 4 He starting from QMC evaluations of the intermediate scattering 
function in imaginary-time |26j; the results were in close agreement with experimental 
data[26j. Moreover the method has allowed to extract also multiphonon energies with 
a good accuracy level. When applied to the two curves depicted in the upper panel 
of Fig. 3.4, we find the two spectral functions in the lower panel of Fig. 3.4| it is 



apparent that the lowest-w peak is indeed well resolved from higher-energy contri- 
butions. Crucially, its position does not depend on the actual choice of the operator 
Af, and it is in excellent agreement with the smallest decay constant found by the 
simple exponential fit. The spectral weight instead is different, consistently with the 
differences between PW-FC and BF-FC. 



3.0.5 Results 



We fit a fifth order polynomial to the density dependence of the energies of the 
triangular crystal and of the paramagnetic and the ferromagnetic fluids, listed in 
Table 3.0.5 The resulting equation of state of two-dimensional 3 He is shown in 
Figure 3.0.5 With the fermionic correlations method, we find a transition between 
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Table 3.4: The equations of state of 3 He for the paramagnetic fluid and the solid 
(solid lines in Figure 3.0.5) are of the form a±p + a 2 p 2 + ct3p 3 + «4p 4 + c^p 5 . This 
Table lists the values of the parameters Oj, with lengths in A. 





liquid 


solid 


ax 


21.23782 


57.35474 




-1344.413 


-2598.784 


a 3 


45093.37 


58695.29 


OLA 


-569306.0 


-532201.7 


«5 


4383507 


3063129 



4 




Figure 3.5: Equation of state of 3 He in two dimensions. Solid line (broken across the 
coexistence region): liquid and solid 3 He; dashed line: mass-3 boson fluid; dotted 
line: liquid 3 He, after Ref. [12]. The latter is only reliable at low densities. 



the paramagnetic fluid and the triangular crystal around p = 0.061 A -2 , with a 
narrow coexistence of about 0.002 A~ 2 , while the ferromagnetic fluid is never stable 
(see Table 3.5). The obtained solidification density can not be directly compared 
with experimental data since we are studying a model of an ideal 2<i-system. It 
could be interesting in future calculations to consider an adsorbing external potential 
representing the interaction of the 3 He atoms with a substrate. 

The energy of the bosonic mass-3 liquid is also shown. This fictitious system, 
simulated to extract the PW-FC and BF-FC energies, crystallizes at p = 0.069 A~ 2 . 
The freezing density of 3 He is considerably higher than the highest density simulated 
in Ref. [12]. Correspondingly, the equation of state given in Ref. [12] is only reliable 
at relatively low density. In particular, while it is only slightly below our results for 
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Table 3.5: Ground state energy of 3 He in K, calculated by the FC method for the 
fluid phases and assumed to equal the bosonic energy for the solid phase. 



density 


liquid C = 


liquid C = 1 


solid 


0.020 


n 1 7A7/ 1 o\ 

0.1707(18) 


0.3218(25) 




U.U45 


0.8168(86) 


0.9075(86) 




0.050 


1.1500(ol J 


i.21zo(oo ) 




0.055 


1.5972(93) 


1.6574(91) 




O.OoO 


2.2Uoy(oo J 


2.24yo(54) 


2.2o0b(o4) 


O.Ooo 


O HAS f TO \ 

o.OOoo(fo) 


O.0o5y(4o ) 


2.yiyo(2o ) 


0.070 


4.0644(33) 


4.0915(34) 


3.7878(35) 


0.075 






4.8728(44) 


0.080 






6.2445(35) 


0.085 






7.9589(39) 


0.090 






10.0661(46) 


0.095 






12.6739(39) 


0.100 






15.8536(45) 



p < 0.045 A -2 as a consequence of the difference of interparticle potential adopted[ 
it becomes even lower than the bosonic equation of state near the melting density, by 
an amount far larger than what could be due to the different employed potential. 

The treatment of the spin polarization state requires a special care [23 EJ ESI 
[9]. In contrast to Ref. [12], we find that the BF-FN energy can be significantly 
higher than the unbiased Fermionic correlations (FC) energy. Starting from negligible 
values at low density, the BF-FN error quickly increases approaching the strongly 
correlated regime. As expected[8], it is larger for the paramagnetic than for the 
ferromagnetic fluid. This happens, because the available variational wave function 
for ferromagnetic states are more accurate than those for paramagnetic states. These 
findings are exemplified in Figure 3.6 The inadequacy of the BF-FN is striking in 
the phase diagram: Figure 3.7 shows that BF-FN incorrectly predicts a transition to 
a ferromagnetic fluid well before crystallization takes place. Such a transition is also 
evident from Figure 3J3, which shows the BF-FN results for the polarization energy 
e(C) at various densities. The unbiased results, shown in Figure 3.9, display instead 
a paramagnetic behavior even in a metastable fluid phase well beyond the freezing 
density. 

From the FC polarization energy e(() we can estimate the spin susceptibility 
enhancement x/Xo- Assuming a quadratic dependence over the whole polarization 
range, which is generally consistent with the data of Figure 3J3, we find an excellent 
agreement with the measured susceptibility. Figure 3.10 shows the comparison be- 
tween the calculated x/Xo an d the experimental data. We display only the results 
obtained in the second layer of 3 He on graphite [3] since they extend to the highest 
density in the fluid phase, but experiments carried on with differently preplated sub- 
strates lead to equivalent results in their respective density ranges. The agreement 
among the results obtained using different substrates induces us to expect that our 
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Figure 3.6: (Color online) Upper panel: Bose-Fermi gap, Abf, as a function of the 
spin polarization, (, at density p =0.020 A~ 2 evaluated via PW-FN, BF-FN, and 
BF-FC with N = 18 particles. Middle panel: Bose-Fermi gap, A BF , as a function of 
the spin polarization, (, at density p =0.045 A -2 evaluated via PW-FN, BF-FN, and 
BF-FC with N = 18 particles. Lower panel: Bose-Fermi gap, Abf, as a function of 
the spin polarization, (, at density p =0.070 A -2 evaluated via PW-FN, BF-FN, and 
BF-FC with N = 26 particles. The statistical uncertainties are below the symbols 
size. 



ideal model actually captures the physical mechanisms underlying the behavior of 
X/Xo- 

The results for x/Xo evaluated with BF-FN calculations diverge at a density 
around 0.050 A -2 , consistently with the BF-FN prediction of a phase transition taking 
place around the above mentioned density. The need for an exact QMC approach 
is thus witnessed by the failure of the BF-FN approximation to predict the lack 
of a polarization transition experimentally observed in the fluid phase, let alone an 
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Figure 3.7: Unbiased FC versus Fixed-Node equation of state. Thick solid line (bro- 
ken across the coexistence region): paramagnetic liquid and solid 3 He (FC); dashed 
line: paramagnetic liquid (FN); dotted line: ferromagnetic liquid (FN); the dashed 
and dotted lines terminate at the FN freezing density; thin solid line: energy of the 
solid, down to the FN melting density. For each density, the energy is relative to the 
energy Eq of the mass-3 boson fluid. 



accurate value for the spin susceptibility. We emphasize that in principle it is possible 
to improve the fixed node results working on the choice of the trial wave function. 
Our purpose in this work was to follow a methodology which is unbiased, that is 
independent on the choice of the wave function; such methodology gives access only 
to the energy and its derivatives. Here we have shown the improvements with respect 
to the results obtained with a particular fixed-node approximation that has already 
been usedp2] for 3 He. 

In conclusion, we have calculated the equation of state and the polarization energy 
of 3 He in two dimensions by means of an unbiased QMC method. The system crys- 
tallizes into a triangular lattice from the paramagnetic fluid at a density of 0.061 A~ 2 , 
with a narrow coexistence region of about 0.002 A~ 2 ; the ferromagnetic fluid is never 
stable. From the polarization energy we obtain a spin susceptibility enhancement in 
excellent agreement with the experimental values. 

We remark that, although the Fermionic correlation technique is in principle un- 
biased, to obtain the estimation of the Bose-Fermi gap one has to face an ill-posed 
inverse problem: the inversion of the Laplace transform at the presence of a limited 
set of noisy data; the quality of the results of the inversion procedure cannot be guar- 
anteed a priori, in this work we have found empirically that the obtained correlation 
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Figure 3.8: (Color online) Fixed-node results for the polarization energy E(Q — 2^(0) 
relative to the Fermi energy E F at p = 0.020 (open triangles), 0.045 (open squares), 
0.050 (filled squares), 0.055 (open diamonds), 0.060 (filled diamonds), 0.065 (open 
circles), 0.070 (filled circles) A -2 , i.e. from top to bottom. The function E &t (() 
is a quadratic polynomial in ( 2 fitted to the simulation data; the solid line is the 
density-independent result for non-interacting particles. 



functions could be safely inverted obtaining robust results, which has been checked 
using different techniques. 

Moreover, the estimation of the Bose-Fermi gap via the Fermionic correlation 
method is limited to relatively small systems: the present results are obtained with 
either 18 or (in most cases) 26 particles. While the size effect remains the main 
source of uncertainty of the present calculation, the agreement of the calculated and 
measured spin susceptibility suggests that finite-size errors are relatively small. 
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Figure 3.9: (Color online) Exact results for the polarization energy E(Q — E^ t (0) 
relative to the Fermi energy Ep at p = 0.020 (open triangles), 0.045 (open squares), 
0.050 (filled squares), 0.055 (open diamonds), 0.060 (filled diamonds), 0.065 (open 
circles), 0.070 (filled circles) A -2 , in order of decreasing dispersion. The function 
-Efit(C) is a quadratic polynomial in ( 2 fitted to the simulation data; the solid line is 
the density-independent result for non-interacting particles. 
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Figure 3.10: Enhancement of the spin susceptibility as a function of the density: 
(filled circles) as measured in the second layer of 3 He on graphite [3]; (open circles) 
as calculated assuming a quadratic dispersion over the whole polarization range in 



Fig. |3.9[ The corresponding Fixed-node result from Fig. |3.8| would diverge at p 
0.050 A" 2 . 
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Chapter 4 

• • • Q 

Dynamics of two— dimensional He 



In the previous chapter we have studied static properties of 2d 3 He; in particular 
we have inspected the dependence of the energy versus spin polarization at different 
densities, showing that 2d 3 He remains a paramagnetic fluid up to the freezing density. 
In this chapter we focus on the dynamical properties of this system. In fact, the low 
energy dynamics of 3 He is of outstanding importance in condensed matter physics to 
understand the thermodynamic behavior of quantum strongly correlated systems pp. 

Recently inelastic neutron scattering experiments have been performed on a mono- 
layer of liquid 3 He adsorbed on suitably preplated graphite: for the first time the col- 
lective zero-sound mode has been detected as a well defined excitation crossing and 
possibly reemerging from the particle-hole continuum typical of a Fermi fluid [TU E]. 
In this chapter, we undertake an ab-initio study of the low-energy collective excita- 
tions, in particular the zero-sound mode, of a strictly two-dimensional (2d) 3 He sam- 
ple relying on Quantum Monte Carlo (QMC) methods. This is particularly appealing 
since it has been shown that the strictly 2d model is often a realistic representation 
of the adsorbed liquid layer, as far as the liquid phase properties are concerned 0, E]. 
The key quantity to be computed to compare with neutron scattering experiments 
is the dynamical structure factor S(q, u), which, apart from kinematical factors, is 
related to the differential cross section and contains informations about low-energy 
dynamics of the sample; in the case of 3 He, the dynamic structure factor is a sum of 
a coherent term S c (q,u) and an incoherent contribution due to the coupling of the 
nuclear spin with the neutron beamjl], Si(q,cu) 

S(q,u) = S c (q,u) + (ai/a c )Si(q,uj) (4.1) 

1 r+oo 

SMoo) = ^- b j dte**{e* B p-j) (4.2) 

i r+oo 

Si&v) = ^J dte**(e i i B p, 9 e- i l B p,_ 9 ) (4.3) 

The brackets indicate a Ground state or thermal average, H is the Hamiltonian 
operator, p$ = J2iLi e~^' n and p z * = J2f=i e~ l ^ r * — $^i=i e~ l ^' r i are respectively the 
local particle and spin densities in Fourier space. The parameter b is the coherent 
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scattering length and a c and Oi are the scattering cross sections for the coherent and 
incoherent scattering. Similarly to the previous chapter, we are interested here only 
in zero temperature properties. The excitations of the system manifest themselves 
in the shape of S(q,u), appearing either as sharp peaks if they are long-lived or as 
broad structures if strong damping is present. In particular, the zero sound mode, 
which is the main target of this work, is related with S c (q,u); the ratio o~i/o~ c has 
been shown[5] to be 0.20(5); moreover, in the experimental data in Ref. [6] there 
is a well defined signal from the zero-sound mode but possible excitations from the 
incoherent part of S(q,u) (i.e. spin waves) are much harder to discern. Given these 
considerations, the data in Ref. [6] is dominated by S c (q,u); we focus here on the 
coherent dynamical structure factor, S c (q,co). QMC methods may indeed give access 
indirectly to the dynamic structure factor, S c (q,oj), because they allow to evaluate 
the intermediate scattering function: 

F%t) = {e rk pie^* p-t) (4.4) 

by simulating a stochastic dynamics in imaginary time driven by the simple Hamil- 
tonian: 



fc2 N N 

H = -^Y.^ + ll v ^-r J ) . (4.5) 

3 i=\ i<j=l 

Here is the mass of 3 He atoms and the pair interaction v(r) is a realistic effective 
potential among 3 He atoms 0. 



Since the ground state is not known, a QMC calculation of (4.4) requires an 
additional time f to project a trial wave function ijjj, onto the exact ground state "0(f 
(see chapter [2]): 



The correlation function (4.4) is the Laplace transform of S c (q,ui). Despite the 
well known difficulties related to the inversion of the Laplace transform in ill-posed 
conditions, the evaluation of S c (q,u) starting from the QMC estimation of F(q,r) 



(4.6) has been proved to be fruitful for several bosonic systems[Bl El ITU] . 

For a Fermi liquid, the difficulty is further enhanced by the famous sign problem, 
thereby the computational effort grows exponentially with the projection time (as 
well as with the number of particles). The total projection time 2f + r in Eq. ( |4.6 ) 



is too large for all practical purposes. 

While accurate approximations exist to circumvent this problem in the calculation 
of static ground-state properties [TT], we are aware of no applications of approximate 
schemes such as the restricted path [12] or constrained path [13] methods to the cal- 
culation of imaginary-time correlation functions. 

We thus resort to the following approximation: 

^ = e~ ffl ^ ~ Ve~ f{i ^ = Vi[)* (4.7) 
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where a superscript F(B) indicates Fermi(Bose) statistics and V is a Slater determi- 
nant. In the resulting approximate correlation function 



F A (q,r) 



(4- 



the projection time between the determinants, which determines the severity of the 
sign problem, is limited to r; Fa is an approximation of the intermediate scattering 



function in imaginary time (4.4), which would be exact if T>ip B were the exact Fermi 
ground state. Its inverse Laplace transform is an approximation of the dynamical 



structure factor (4.2). For a given wave vector, the positions of the peaks in S(q,u) 



provide the energy of the excitations, while their shape is related to the life-time of 
the excited states. In general, the approximation (4.7) introduces biases both in the 



excitation energies and in the shape of S c (q,u). In order to enhance the robustness 
of this approach, we introduce also another correlation function, Fb, which is defined 
on the bosonic ground state: 



F B (q,r) 



B\„tH jy* 



(4.9) 



Despite Fb is not directly related to the dynamical structure factor of the Fermi 
liquid, this function has some useful features: on one hand, it contains the exact 
fermionic spectrum, as can be seen from the spectral resolution: 



+oo 



F B (q,r) = J2' 



-t(E* 



l(p- g -p ^\< )\ 2 

<^w> 



(4.10) 



On the other hand, it is a bosonic correlation function and thus it can be evaluated 
with great accuracy by means of exact bosonic QMC methods. If, moreover, the 



approximation (4.7) is accurate enough, the coefficients b n become, apart from an 



unessential normalization, the spectral weights /„ of the exact intermediate scattering 

\(P-^\0\ 2 



function (4.4) 



fn 



(4.11) 



We note finally that Fb arises as a natural generalization of the Fermionic correlations 
method: in fact, Fb has the same form of Eq. (3.11) that has been used in the 



previous chapter about the ground state of an 3 He film: in that context, the Fermionic 
correlations method provided results for the magnetic properties of the system in 
impressive agreement with experimental data. 

We argue that a comparison between dynamical properties evaluated with Fa and 
Fb might provide a strong indication of the robustness of our approach. 

We studied a system of N = 26 structureless 1/2-spin fermions of mass m 3 , in- 
teracting via the Aziz potential described by ref. [7], a very accurate model for the 
effective interactions among 3 He atoms. The choice of the particle number was in- 
spired by our previous work in Ref. [3]: such number of atoms was chosen to be a 
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closed-shell number; this choice minimizes size effects related to the discrete Fermi 
sea, but still allows to extract physical information from the imaginary time correla- 
tion functions, which rapidly become steeper as the size of the system is increased. 

Differently from the work described in the previous chapter, we have not used 
twisted boundary conditions (TBC). This choice is motivated by the fact that the 
effect of TBC enters in the estimation of both the Fermi-Bose gap and the energy of 
the excited state with respect to the bosonic ground state; being the energy of the 
excitation a difference between these these two quantities, we assumed that, as a first 
approximation, the effects of TBC cancel out. In conclusion, considering that the 
evaluation of the necessary correlation functions are required with high quality data, 
neglecting TBC is a good compromise between accuracy and practical computing 
times. 

We have focused on a density around 0.047 A" 2 , close to the experimental 
conditions |14j. Moreover, we have explored the behavior of the sample at the densi- 
ties 0.038 and 0.060 A -2 in order to investigate the possible density-dependence of 
the excitations of the system. In particular, the highest density was chosen very close 
to the freezing point. The QMC evaluation of Fb requires a simple generalization of 
the methodology that we have followed in the previous Chapter: a fictitious system 
of bosons of mass m 3 is simulated with the Shadow Path Integral Ground State 
method. The imaginary-time propagation was 1.3125 K -1 and the density matrix 
approximation was a Pair Productpj)] with imaginary-time-step of 1/160 K _1 . 

The Shadow Path Integral Ground State [16] technique was chosen in the compu- 
tation of both the bosonic ground-state energy and the correlation functions. 

Performing such a QMC simulation, we have computed Fs{q, t) for each wave- 
vector q, together with the correlation function: 



*b(r) = v,' , /; r o/ (4.12) 



(^\V*e- TH V\^) 

WW) 

which is precisely the correlation function that was used in the previous chapter 
in order to estimate the energy gap between the bosonic fictitious system and the 
fermionic ground state. Fa has been then estimated from the exact identity: 

FMr) = F 0^ (4.13) 

It is well known that, in order to extract information from imaginary time cor- 
relation function, an inversion of the Laplace transform in ill-posed conditions is 
necessary. This can be carried out by means of the Genetic Inversion via Falsification 
of Theories (GIFT) [5], which has already provided very accurate results in the study 
of low energy excitations of Bose superfluidsp] and supersoridsfTO]. Naturally the 
problem is unavoidably ill-posed: the quality of the results of the inversion procedure 
cannot be guaranteed a priori; however, a test of reliability of the inversion procedure 
can be obtained by comparing our estimations of the dynamic structure factor with 



experimental data; Fig. 4.3 shows a remarkable agreement. 
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Figure 4.1: Comparison between the spectral functions obtained from Fa and those 
obtained from Fb for some wave-vectors q. The two spectral functions have a com- 
patible shape, with a shift in energy of Eq — E B . The data for Fa-, differently from 



Fig. 4.3, has been obtained from the average of six different evaluations of Fa(t) 



In Fig ]4.1| we show a comparison between the estimated inverse Laplace transforms 
of Fa and Fb ■ apart from a rigid shift in energy, due to the difference Eq — E B which 
we have estimated in the previous chapter, the reconstructions coincides within the 
"algorithmic resolution" of the GIFT methodology. This represents a confirmation 
for the robustness of our approach. We remark that Fb is much easier to handle than 
Fa, since it does not suffer of long-r large fluctuations due to the presence of the 



r-dependent denominator in (4.13). 



The shape of the reconstructed spectral functions depends on the lowest-energy 



fermionic exact eigenstates not orthogonal to the wave function f)-^T> ip B , a state 
containing a density modulation of wave vector q. 



Our assumption (4.7) asserts that such wave function is very similar to p^fipQ. 



As appears evident in Fig, 4. 2 at low wave vectors the inversions of Fs(q, r) provide 



spectral functions with sharp peaks; this provides our microscopic estimation of the 
zero-sound mode dispersion relation, with an "algorithmic resolution" similar to that 
found in the previous chapter, where the Fermi Ground State signal was detected, at 
higher wave vectors the peaks become much broader. We interpret such a broadening 
as a damping of density fluctuations due to the presence of other excitations, in 
particular the particle-hole excitations. 
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Figure 4.2: Panel a) the gap between the B osonic and the Fermionic ground state 
(dashed line) and the inversion of eq. ( |4.10[ ) for q = 0.0267 A" 1 (filled line). Panel 
b) the same of panel a), but with q = 0.801 A -1 



In Fig,4 1 3 we show the comparison between our estimation of the dynamic struc- 
ture factor of the 3 He film and the experimental datajH]. The agreement is impressive 



and gives a strong support to the approximation (4.7); it is also clear from this com- 



parison that Eq. (4.7) describes accurately also the mechanisms which give rise to a 



broadening of the dynamic structure factor; this is displayed even better in Fig, 4.4 



where we report, in a color plot, the estimated Sb(<1, oS). At low q we find well defined 
excitation energies, while, as the wave vector increases, we observe the sharp mode 
becoming damped. 

What kind of excitation provide the damping of the zero-sound mode? Any expert 
in Fermi liquid theory would immediately answer: the particle-hole continuum. But 
what does it mean in an ab-initio approach to a strongly correlated fermion fluid? 
Our idea is to exploit the correlation functions formalism to build up a "Fermi-liquid 
like" function: 



ph\ 



* p~ tH T) l. 



(4.14) 

where T> p h is simply a Slater Determinant like that employed in the previous 



chapter (see Eq. (3.15)), with Q\ = 0. Differently from Eq. (3.15), in the enumeration 



of the wave-vectors {k n } } one element has been taken out of the ideal gas Fermi sea; 
the bosonic ground state and the backflow correlations provide a "dressing" for the ph- 



wave function T> p hipQ . In Fig. 4.5 we show both the particle-hole excitation energies 



for the ideal Fermi gas, which do not form a continuum since the system is finite, 
and the estimated energies extracted from the inversions of F p h(r), which have been 
evaluated for wave vectors at the high-g borderline of the structure; we focused on 
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Figure 4.3: From left to right the obtained dynamic structure factor for increasing 
wave vectors at p = 0.047 A~ 2 . The yellow shadow represents statistical uncertainties 
obtained from six different evaluations of the dynamic structure factor at each wave- 
vector; filled circles are the available experimental data from Ref. ([II]). The wave- 
vector shown in picture are those accessible from our simulation, the experimental 
wave vectors are q = 0.55 A -1 (b), q = 1.15 A -1 (d), q = 1.25 A -1 (e) and q = 
1.65 A- 1 (f). 



the high-g borderline in order to verify whether the roton states reemerges from the 
particle-hole band or not. The particle-hole band of the ideal Fermi gas has been 
computed following the definition of particle-hole energy: from an ideal Fermi gas of 
N particles of mass m in a square box of late L, the Fermi Sea is labeled by quantum 
numbers that define the wave vector of each state, ky = (2tti/L, 2irj/L); a particle- 
hole excitation characterized by a hole at (io,jo), corresponding to a wave-vector k Q 
inside the Fermi sphere, and a particle at (ix,ji) with wave-vector k\ outside the 
Fermi sphere, has a wave-vector q and an energy E l ^ eal defined as: 



E 



ideal 



(E F + 



2m 



(Ep — 





2tt 


Q = 




h 2 




2 




ki 


) 


2m 





V(*i _ ^o) 2 + O'l - jo) 
2m 



ko 



(4.15) 
(4.16) 



where Ep is the Fermi energy of the system. The particle-hole band of the ideal Fermi 
gas in Fig. 4.5 has been obtained following this prescription; all the possible particle- 
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hole combinations which gave a wave vector q in the displayed range were considered. 
Comparing the ideal particle-hole band with that of the interacting system one can 
also estimate the effective mass m*; in this context, the effective mass is a parameter 
of the Landau Fermi liquid theory which gives the mass of the quasi-particle; m* 
is obtained from the ratio of the particle-hole energies of the interacting and the 
non-interacting systems at the same wave-vector, 



E ideal = « |^2 (4 



2m 

E T = kf (4-18) 



pint | 2 

2m* 

j^jideal 



m 



m E 1 ™ 1 



(4.19) 



This relation, however, is valid only in the range of applicability of the Landau Fermi 
liquid theory; in particular, this evaluation of the effective-mass holds for small wave- 
vectors; in this work we give a rough estimate of the effective mass as the average of 
E idealj E int for each part i c l e -hole wave-vector q computed in our simulations. In the 
interacting system, we find in general that the particle-hole energies become smaller, 



resulting in an higher effective mass (see Tab. 4.1). In particular this has important 
consequence as far as the re-emergence of the zero-sound mode from the particle-hole 
band is concerned: in contrast to what the authors of Ref. [B] argue, using the non- 
interacting estimation of the particle-hole band, we do not observe such re-emergence 



in the roton region at any density (see Fig, 4. 2) 



Table 4.1: Effective to bare mass ratio estimated from the computed particle-hole 
excitations in a system of N = 26 particles at the studied densities. 



Density ( A~ 2 ) 


m 


0.038 
0.045 
0.060 


1.3(3) 
1.8(1) 
2.0(1) 



We point out that this evaluation of the particle hole excitations gives only a first 
evidence that the roton mode is still inside the particle-hole band; a further step on 
this topic consists in a size scaling analysis that has been planned for future works. 

The static response function can be obtained from the knowledge of the dynamic 
structure function. The described method can thus be viewed also as a new and very 
accurate way to compute the static response function of a fermion system. From 
six independent evaluations of the dynamic structure factor we computed the static 
response function defined as \q — ~ 2p / du S{ ~^ and the result is displayed in figure 



4.3 To our knowledge, this is the first microscopic ab-initio computation of the static 
response function of two-dimensional 3 He. 

As a last result, we note that Eq. (4.8), for r = is the definition of the static 
structure factor evaluated on a Fermi state [D^q). This state, as shown in Fig. 4.3 
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is a good approximation of the true fermionic ground state; this is at least true for 
the low-energy dynamical properties; we assume that this holds also for the static 
structure factor. In Fig. 4.6 we show the static structure factor of 2d 3 He compared 
with that of the fictitious "bosonic" 3 He. The similarity between the two static 
structure factors is evident; this indicates that the structural properties of Helium 
are dominated by the inter-atomic potential rather than the quantum symmetry. 

In conclusion, we have presented the first ab-initio computation of the zero-sound 
excitation energy. We have also proposed an approximate evaluation from first prin- 
ciples of the dynamic structure factor that is found to be in very good agreement with 
experimental data|14] . We employed a well tested methodology involving the Laplace 
inversion of imaginary-time correlation functions [TTl [3] and extended it in order to 
handle the excited states. Our results are in agreement with the experimental data 
and show that our variational estimation of the dynamic structure factor is accurate 
enough to represent the broadening of the zero-sound mode in the particle-hole band. 
At the studied densities we did not observe the re-emergence of the roton mode from 
the particle-hole band, however our data on particle-hole excitations is not yet con- 
clusive: possible finite size effects on the particle-hole band have still to be studied 
with simulations of bigger systems. 
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Figure 4.4: The horizontal continuous line represents the energy gap between the 
Fermionic and the Bosonic ground state, slightly above - at wave vectors ranging 
roughly between 1 and 2 A -1 - the discrete particle-hole right boundary is shown; 
the remaining and larger vertical bands are the density-density collective excitations. 
The bands are centered on the corresponding values and their width has been enlarged 
for a better visibility. 
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Figure 4.5: (Circles) Particle-hole excitations for N=26 non-interacting atoms of 
4 He mass. (Dashed lines) Particle-hole band for the ideal gas in the thermodynamic 
limit. (Filled circles) Particle-hole excitations for N=26 atoms of interacting 4 He. 
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Figure 4.6: The static response function of 3 He obtained from eq. (4.10). (Circles) 
p = 0.038 A" 2 . (Triangles) p = 0.047 A~ 2 . (Squares) p = 0.060 A~ 2 . 
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Figure 4.7: The Static Structure Factor of 2c? 3 He (filled symbols) compared with that 
of "bosonic" 2d 3 He (empty symbols). Data is relative to three densities: p = 0.038 
A" 2 (Circles); p = 0.047 A" 2 (Squares) and p = 0.060 A" 2 (Triangles). The dashed 
lines are guides to the eye. 
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Chapter 5 



Study of He adsorbed on 
Graphene— Fluoryde and 
Graphene— Hydrate 

At the forefront of current research in condensed matter physics is the study of 
strongly interacting systems, with a remarkable variety of phase transitions pp. The 
effects of fluctuations are enhanced in low dimensions and in the presence of frustra- 
tion [2] . These represent some of the motivations for studying adsorption phenomena, 
where important roles are played by the gas-gas interaction and the "tunable" effect 
of the substrate. The surface of graphite has long been a playground for studying 
two-dimensional (2D) monolayer phases of classical and quantum gases [3]. 

Probably the best understood adsorption system is the He monolayer on graphite 
[I]. Experiments carried out at the University of Washington ca. 1970 revealed for 
the first time behavior corresponding to a two-dimensional (2D) gas. More dramatic 
was the appearance of a spectacular peak in the specific heat of 4 He near T c = 3 
K. This peak, well described by the 3 state Potts model, manifested a 2D transi- 
tion from a high T fluid to a low T commensurate (a/3 x a/3 R30°) phase, provid- 
ing a benchmark measure of coverage, not seen in previous adsorption experiments. 
This ordered phase (at density = 0.0636A~ 2 ) corresponds to atoms localized on 
second-nearest neighbor hexagons. At higher densities near completion of the first 
monolayer (p = O.llA -2 ) an incommensurate 2D triangular solid phase is present; the 
phase diagram at intermediate densities is not yet completely determined. A quan- 
titative understanding of the He-graphite interaction was made possible by precise 
scattering measurements of surface bound states and band structures [S1E]- 

The availability of graphene (Gr) and its derivatives like graphane (GH) [7] and 
graphene-fluoride (GF) j8] offers the prospect of novel adsorption phenomena. 

Since Gr is just a single plane of graphite, the symmetry and corrugation are 
expected to be very similar in the two cases. If Gr is rigid, no new phenomena are 
expected for adsorption on one side of Gr, in comparison with graphite, [5] and this has 
been verified by recent quantum simulations of 4 He [8]. The situation is different for 
the derivatives of Gr, graphene-fluoride (GF) [7j and graphane (GH) [TJ |9] that have 
been recently obtained experimentally. Because GF and GH have surface symmetries 
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and compositions which are quite different from Gr, adsorbed gases will have very 
different properties. 

In the next section we show a model adsorption potential for He on GF and GH. 
We will then show the study of a single 4 He and 3 He atom on these substrates, as 
well as submonolayer films of 4 He at coverages similar to that (p = 0.064A -2 ) of the 

x \/3 R30° state on graphite. In the last section the properties at high coverages 
will be described. 



5.0.6 Adsorption potential 

Graphane and graphene-fluoride have a similar geometry; half of the H (F) atoms 
are attached on one side of the graphene sheet to the carbon atoms forming one of 
the two sublattices of graphene. The other half are attached on the other side to the 
C atoms forming the other sublattice. The H (F) atoms are located on two planes 



see left of Fig. 5.2); one is an overlayer located at a distance h above the pristine 



graphene plane while the other is an underlayer at a distance h below the graphene 
plane. In addition, as seen in Fig. |5.1| there is a buckling of the C-plane with the 



C atoms of one sublattice moving upward by a distance b while the other sublattice 
moves downward by the same amount. A He atom approaching GH (GF) from above 
will interact primarily with the H (F) overlayer, but it will interact also with the C 
atoms and the H (F) atoms of the underlayer. 




=C ~)=F(H) 



Figure 5.1: Geometry of the substrate with the definitions of the buckling parameter 
b, the interplane distance h and d, the carbon-carbon distance on the plane 



We have adopted a traditional, semi-empirical model to construct the potential 
energy V(r) of a single He atom at position r near a surface PUJ HH [121 US]- The 
potential is written V(r) = V iep (r)+V ai ,t{'^), a sum of a Hartree-Fock repulsion derived 
from effective medium theory, and an attraction, V a tt(f), which is a sum of damped 
He atom van der Waals (VDW) interactions and the polarization interaction with 
the surface electric field. The first term is V rep (r) = ap(r). Here a = 364 eV-bohr 3 
is a value derived by several workers as the coefficient of proportionality between 
the repulsive interaction and the substrate's electronic charge density p(f) prior to 
adsorption. The geometry of GH and GF, their electronic charge density and the 
electrostatic potential have been obtained using Density Functional Theory with an 
all-electron triple numerical plus polarization basis set with an orbital cutoff of 3.7 A 
as implemented in the DMol3 code[14j. The exchange and correlation potential was 
treated in a Generalized Gradient Approximation parameterized by Perdew, Burke, 
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Figure 5.2: Two schematic views of GF. F (C) atoms are light (dark) gray. Positions 
of atoms are to scale but their sizes are arbitrary. The black balls represent two 
adsorption sites for He, one of each kind. GH is similar. 



and Ernzerhof [T5|. We use a tetragonal unit cell containing four C atoms and four 
H (F) atoms for GH and GF, respectively. The cell dimensions for GF are at = 2.59 
A, ci2 = 4.48 A, and 03 = 12 A, while for GH we use at = 2.52 A, ci2 = 4.36 A, and 
03 = 12 A. The Brillouin zone was sampled with a Monkhorst-Pack grid of 6 x 3 x 1 
k points in both cases. The self-consistent cycles were run until the energy difference 
was less than 10~ 6 eV. The atomic positions were relaxed until the forces on all atoms 
were lower than 0.01 eV/A. As a result, the C-F distance is 1.38 A, the C-C distance 
1.57 A, the C-C distance projected on the x—y plane is d = 1.495 A and the buckling 
displacement b = 0.484 A; while in GH, the C-H distance is 1.11 A, the C-C distance 
1.52 A, d = 1.453 A and b = 0.45 A. 

The attraction is a sum of contributions; for GH, 

Ktt(r) = V n+ (r) + V gr (f) + V H _(f) - «HeE 2 (f)/2 (5.1) 

The right-most term is the induced dipole energy, where an e = 0.205 A 3 is the 
static polarizability of the He atom and E(r) is the electric field due to the substrate. 



This term gives a minor contribution to the adsorption potential (see Fig. 5.3) and 
has been neglected in this work. 

The three VDW terms for GH originate from the H overlayer, the graphene sheet 
(we are neglecting in this term the small buckling of the graphene sheet) and the 
H underlayer, respectively. These terms are described by the attractive part of a 
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Figure 5.3: Upper figures: 4 He on GF. (Left) minimum value with respect to z 
of V a d s {f) along the direction (x,0) (black) and (0,y) (red). Full lines represent the 



adsorption potential, V a , obtained using Eq. (5.1 ) for the attractive part, dashed lines 
the adsorption potential, V&, obtained neglecting the induced dipole energy. (Right) 
the relative difference (in percentage) between V a and Vb with respect to V a , namely 
100 * (Vb — V a )/V a . Lower figures: the same for GH. 



Lennard- Jones potential, 



V H -(r) 
V gr {r) 



-E 



3 \ r ~ r j 



H- 



-E 

j 

E 



6C 



a 



6h 



J \ r ~ r 3 



(5.2) 
(5.3) 
(5.4) 



where the sum spans over the carbon or hydrogen positions; Cqc an d Cqh are respec- 
tively the Helium-Carbon and the Helium-Hydrogen VDW coefficients. Due to the 



distance between the helium monolayer and the graphene plane, Eq. (5.2) can be 
integrated over the x-y plane 



-6h-Cqh 



d 2 R- 



1 



z + h) 2 + 



R 



A 



c 



(z + hy 



(5.5) 
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where Oh- is the Hydrogen density of the sublayer on the x-y plane and h = h + b/2 
is the distance between the Hydrogen underlayer and the mean Carbon plane (see 



Fig. 5.1). The He-H VDW coefficient, Cqh, is obtained from the VDW interaction of 
He-H 2 . In Ref. [IE], the potential energy between He and H 2 is written as a sum of 
potential terms regarding the interaction of Helium with each single Hydrogen, 

U He - H2 (r, R u R 2 ) = U He - H {\r- R x \) + U He . H (\f- R 2 \) (5.6) 

where R\ and R 2 are the positions of the Hydrogen atoms. Considering an attractive 
VDW term, —Cqh/^i f° r each Hydrogen atom and neglecting the structure of the 
H 2 molecule, we approximate the attractive part of UHe-H 2 as an isotropic VDW 
interaction C6# 2 /r 6 . We then have: 

Ce H = (5.7) 

The value of Cqh 2 is obtained from Ref. [T7j; in that work the anisotropy of Uue-H 2 
is also studied. The work shows that the leading term of the long range attractive 
part of the He-H 2 potential is 

Ugr-n&i) = ~ E ife t 1 + r *^(«»7)] (5-8) 

n>3 l-^l n 

where R is the distance from the center of mass of H 2 and 7 the angle between R 



and the axis of the molecule. The leading term of Eq. (5.8) is exactly the VDW 
term C 6 ^ 2 /r 6 ; the anisotropy is described by the Legendre term P 2 ; in particular, 
Ref. [IB] shows that the leading anisotropic correction involves T 6 which is of order 
0.1. This quantity, even though might have some relevance, has been neglected: its 
effects can be taken into account; however, the results in this chapter may not change 
qualitatively; this has been checked with an arbitrary change of the VDW parameters 
for the interaction potential of both He-GF and He-GH. Moreover, this work is based 
on a semi-empirical adsorption potential and it is not aimed to obtain quantitative 
data. 

The term V gi: (f) can be treated in an analogous way of the term Vn-(r), 

V gr (rl = -e 9r C 6C [ *R t 1 A 3 = - (4^) = (5.9) 

z 2 + R 



This approximation, however, neglects the buckling of the Carbon atoms and gives 



rise to variations of 8% in the case of GH (see Fig. 5.4). Apart from a shift in energy 



we don't expect that this approximation would change the qualitative behavior of the 



system; in fact, as can be seen in Fig. 5.4, the larger differences are at the maxima of 



the adsorption potential, where the Helium density is lower. Nevertheless, in this work 



we preferred to use directly Eq. (5.3 ). The Helium-Carbon VDW coefficient appearing 
in this equation, Cqc, can be determined from the known[12j VDW coefficient C3 = 
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Figure 5.4: Upper figures: 4 He on GF. (Left) minimum value with respect to z 
of V a ds{f) along the direction (x,0) (black) and (0, y) (red). Full lines represent 
the adsorption potential, V a , obtained using Eq. ( 5.3[ ), dashed lines the adsorption 
potential, obtained with the approximation (5.9). (Right) the relative difference 
(in percentage) between V a and Vj, with respect to V a , namely 100 * (V& — V a )/V a . 
Lower figures: the same for GH. 



180 meV-A 3 for an Helium atom interacting with a half-space of graphite, V(z) ~ 
— C%z~ 3 . This is connected to the VDW interaction ~ — C^c r & by an integral over 



the half-space; in a way similar to Eq. (5.5) we find the relationship between C3 and 



7T 

6" 



nC, 



6C 



(5.10) 



where d = 3.4 A is the distance between two carbon planes of Graphite; n = 6/d 



is the three dimensional density of the half plane of Graphite. Using Eq. (5.3) and 
Eq. fl5~l0| one obtains that A c = 3C 3 d =1.84 eV-A 4 . 

The term in Eq. (5.4), Vn + (r), gives the main attractive contribution; not only 



this term can not be integrated along the x-y plane but its proximity with the Helium 
monolayer requires the use of damping; we have adopted the Tang-Toennies damping 



procedure for this situation [13] ; the term in Eq. (5.4) thus becomes 



Vu+{r) ^ V damp (\f - fj\) 



Vr, 



damp 



X) 



a 



6H~ 



1 _ „-P* V 6 iMl 
x c Z^n=0 n! 



(5-11) 



The function Vd amp has an asymptotic x 6 behavior; for small values of x the diver- 
gence is cured; this can be seen with a Taylor expansion of e~ l3x . Following Ref. [13], 
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the parameter /3 is the decay coefficient in the asymptotic charge density Ph{t) ~ e~ l3r 
due to the H atom; for Hydrogen = 3.78 A -1 . 

In the case of GF, K, tt (f) has an expression similar to Eq. (5.1) with Vu + and 
Vh_ replaced by Vp+ and by Vp_. The same procedure used for GH applies with the 
coefficient C 6F = 4.2 eV-A 6 as given by Frigo et al [32], /3 = 3.2 A" 1 and A F = 1.1 
eV-A 4 . 



Table 5.1: Parameters for the adsorption potential of He on GH and GF 



Parameter 


Value 


Type 




4.2 eVA 6 


GF 


Cqh 


1.206 eV A 6 


GH 


Cqc 


3.447 el/A 6 


GF/GH 


A F 


1.1 e\/A 4 


GF 


A H 


0.3455 el/ A 4 


GH 


F 


3.2125 A" 1 


GF 




3.77945 A" 1 


GH 


7 


53.9392 eVA 3 


GF/GH 



The obtained adsorption potential relies on the electron density of the substrate. 
This quantity is the output of a DFT computation and is in the form of a 3d table 
formatted as iSx, jSy,kSz, value, where 5x,5y,5z is the spatial discretization of the 
electron density table. The adsorption potential that enters as input in the QMC 
simulations is thus a 3c? table; this table is read in the simplest way: if the simu- 
lation box is a cube in the region with positive coordinates, a position r = (x, y, z) 
corresponds to a bin in the table defined by (a,b,c), where a = x/dx, b = y/dy and 
c = z/dz. 



With such model potentials the adsorption sites (see Fig. 5.6) are above the centers 
of each triplets of H (F) atoms of the overlayer, forming a honeycomb lattice with 
the number of sites equal to the number of C atoms, twice as many as those on Gr. 
Half of the sites are above H (F) of the underlayer but the difference between the 
well depths for the two kinds of adsorption sites is very small, below 1%. For GF the 
well depth is 498 K and for GH it is 195 K. These values do not include the induced 
dipole energy which gives a contribution below 1%. The inter-site energy barrier is 
24 K for GF and 13 K for GH. Both values are significantly smaller than the barrier 
height 41K for graphite. In this last shown in Fig. |5.7[ the energy barrier 



does not depend much on the direction in the x — y plane whereas in the case of GF 
and GH the ratio between maximum and minimum barrier height in the x — y plane 
is of order of 4-5: the energy landscape of the two last substrates is characterized 
by a very large corrugation with narrow channels along which low potential barriers 
are present. The motion of the He atom, especially in the case of GF, essentially 
visits only these channels , as though the atom moves in a multiconnected space; 



this is seen in Fig. 5.0.6 Another significant difference is that the distance between 
two neighboring sites is 1.49 A for GF and 1.45 A for GH whereas it is 2.46 A for 
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graphite and for Gr. Prior to these studies, graphite was believed to be the most 
attractive surface for He, with a well-depth a factor of 10 greater than that on the 
least attractive surface (Cs). The present results reveal GF to replace graphite, since 
its well is a factor of 3 more attractive. 




Figure 5.5: Upper panel: local density p(x,y) (in units of A -2 ) as function of x-y for 
4 He at equilibrium density on GF. Lower panel: local density p(x = 0, y) (in units of 
A~ 2 ) in the unit cell (with side a = 4.486 Ain the GF case and a = 4.347 Ain the 
GH case) along the y direction for 4 He at equilibrium density on GF and GH; note 
the logarithmic scale used for p(x = 0,y). Error bars are below the symbol size; lines 
are guides to the eye. 



5.0.7 The QMC parameters 

The computations in this chapter are based on the Path Integral Ground State [20] 
(PIGS) and the Path Integral Monte Carlo (PIMC) methods^!]. As widely explained 
in the first chapter, with these methods we can compute quantum averages of the sys- 
tem at respectively zero and finite temperature; the PIGS method uses the quantum 
evolution in imaginary-time r of a trial wave function If \l/ t is not orthogonal 
to the ground state, and r is sufficiently long, the quantum evolution purges from 
ty t the contributions of the excited states, yielding the ground state energy and wave 
function. A valuable feature of the PIGS method is that it is exact, in principle; the 
results are independent of ^ t [22] and systematic errors may be reduced below the 
statistical uncertainty. The PIMC method applies the Path Integral formalism to the 
quantum thermal average expressed in coordinates representation; this expression is 
then evaluated with Monte Carlo methods; these methods are extensively explained 
in Chapter [2] 

Both the zero and finite temperature simulations on GF and GH have been per- 
formed with the eight order Multi Product Expansion of the small imaginary-time 
propagator; the imaginary-time discretization is 5r = 1/160 K" 1 , which gives a suf- 
ficiently accurate approximation of the propagator; an example in the case of GF is 



given in Fig. 5.8 Due to the computational complexity of PIMC, especially at low 
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0.5 1 1.5 2 2.5 3 3.5 4 

y (A) 

Figure 5.6: Upper panel: plot of the minimum value with respect to z of the 
adsorption potential He-GF, V(r), in K as function of x-y. Lower panel: the same 
for GH. 



temperatures (i.e. 0.5 K), the simulations at finite temperature performed in order 
to obtain an estimation of the superfluid density were carried out with St = i K _1 ; 
this choice does not lead to quantitative results, however we note once again that 
this work is based on a semi-empirical adsorption potential; moreover, a quantitative 
study of the superfluid fraction would require extensive size scaling in order to de- 
termine the normal-to-superfluid transition; the aim of these simulation is thus to 
give a preliminary estimate of the superfluid fraction in order to determine if there is 
superfluidity rather than the exact density. 
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Figure 5.7: Energy barrier in GF, GH and graphite as atom moves along a line 
making an angle 9 with the horizontal direction and following the height z(x, y) 
giving the minimum of V (f) . Plotted energy is relative to energy at the adsorption 
site. 



The trial wave function ty t that we have used in PIGS is the product of a Jastrow- 
McMillan wave function and a Gaussian along the ^-direction 

^ = e" Ef<J=1 fe) e- A ^=i^"^ 2 (5.12) 

where N is the particle number and = \ fi — fj\ is the distance between two atoms 
labeled i and j. he Jastrow parameters are b = 2.84A and m = 5. The Gaussian along 
the z-direction (i.e. the direction perpendicular to the substrate plane) was used only 
far away from the layer promotion density; its parameters have been obtained with 
a fit of the density along the ^-direction, for GF A = 5.6A -2 , z = 3.72A, for GH 
A = 3.0A- 2 and z = 3.85A. 

At high densities, where the probability to occupy the second layer is not negligi- 
ble, a Jastrow wave function has been used, 




(5.13) 



With these trial wave functions, we allowed a 8r = 0.15 K _1 imaginary-time pro- 
jection before computing the ground-state expectation values. The total imaginary- 
time sampled in our calculations was r =0.4 K -1 . The value of r has been chosen 
following a convergence test of the total energy versus the imaginary-time projection. 
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Figure 5.8: Convergence of the total energy per particle versus time-step St for a 
system of N = 26 atoms of 4 He on GF at equilibrium density The horizontal line 
is the convergence value taken as the average of the energy at the three smallest 
timesteps. The green circle represents the used time-step, St = 1/160 K -1 . 



The Worm algorithm [23] was used at both finite and zero temperature respectively 
for the sampling of the permutations and the computation of the one body density 
matrix. 

The computations required on average 10 5 Monte Carlo steps, the heaviest compu- 
tations were those made for the superfluid fraction at zero temperature and required 
approximately 10 7 Monte Carlo steps. 



5.0.8 A single Helium atom on the substrates 

We computed the exact ground state energy of one 4 He atom or one 3 He atom on 
GF and GH, see Table 5.2 The binding energy on GH is similar to that on graphite, 
whereas that on GF is about three times that on graphite. In both cases the ground 
state is delocalized over the full substrate and both kinds of adsorption sites are 
occupied with comparable probability. 

We have also computed the density-density imaginary time correlation function 
in Fourier space; in the case N — 1, at an imaginary-time r, this function takes the 
form: F(k,r) = (p%(t)p_z(0)), Pk( T ) = ex P ik ■ r{r) . Here r(r) is the position of 

the atom at imaginary time r. F(k, r) contains informations on the excited states of 
the system; these informations can be extracted through an inversion of the Laplace 



82 



25 



20 







— T.B. fit: 


?1 


= 5.695, y 2 


= -0.33 






--- T.B. fit: 




= 3.63, y 2 - 


-0.16 




°v 











>1 

g 10 

LU 





25 

20 



'S--.0 



\ / 
\ / 
A 

% / \ / 

K « 
' v \ 



r 



K' 



M 



b) 



--- T.B. 


fit: Yl 


= 6.434, y 2 - 


-0.04 


--- T.B. 


fit: Yl 


= 4.514, y 2 = 


-0.021 




j. 




Figure 5.9: Panel a) Energy of the first band along two directions of the first Brillouin 
zone for 4 He (triangles) and 3 He (circles) on GF. Data along TK beyond the Dirac 
point K give the results in the II Brillouin zone. The dashed lines are fits made with 
the tight binding model on a honeycomb lattice with nearest neighbor parameter 71 
and next nearest neighbor parameter 72 as in legend. Panel b) Same as for panel a) 
for 4 He and 3 He on GH. 



transform that gives the dynamic structure factor S(k,u): 

F(k,r) = j ' due- UT S(k,u). (5.14) 

However, F(k, r) is known only at discrete imaginary-times r m with a statistical 
uncertainty; the inversion of the Laplace transform in such conditions is an ill-posed 
inverse problem; as consequence, the quality of the extracted informations can not 
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Table 5.2: Kinetic, potential and total energies for the ground state of He on GF, 
on GH and on graphite. In the last column the bandwidth A is shown. Numbers in 
parentheses represent statistical uncertainty in the last digit. 



System 


E ki n (K) 


E pot (K) 


Etot (K) 


A(K) 


4 He+GF 
3 He+GF 


46.78(4) 
51.08(1) 


-422.94(1) 
-413.41(1) 


-376.15(2) 
-362.33(1) 


9.6(1) 
13.7(1) 


4 He+GH 
3 He+GH 


20.51(1) 
22.53(2) 


-153.58(1) 
-149.50(1) 


-133.06(1) 
-126.97(2) 


13.6(4) 
19.4(4) 


4 He+Gr 
3 He+Gr 


25.30(4) 
27.05(2) 


-168.49(1) 
-162.87(1) 


-143.19(4) 
-135.82(2) 


9.6(2) 
15.7(4) 



be guaranteed. The inversion of the Laplace transform has been computed with 
the GIFT method explained in Ref . [23] . Basically, the GIFT method uses a Genetic 
Algorithm to explore a space of solutions {S n (k, w)}; the solutions that can reproduce 
F(k, t) with an user-defined accuracy are averaged together to give the solution. In 
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Figure 5.10: Left side: density-density correlation functions F(k, r) for a wave vector 
k corresponding to the point M in the first Brillouin zone (see Fig. 5.9). On the right 



side the respective inversions of Laplace transform, S(k,u>), are shown. 



Fig. 5.10 an example of F(k,r) and its S(k,u) obtained with GIFT is given: in 
these functions F(k,r), the main contribution comes from the lowest energy band; 
moreover, the excitation appears in S(k, cu) as a well defined peak; it is thus possible 
to obtain the energy spectrum of the first energy band; we interpret the width of 
the peaks as the uncertainty associated to the excitation energy at that wave-vector. 
The computed energy spectrum along the directions TK and TM for He on GF and 
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on GH is shown in Fig. 5J) These bands are represented rather accurately by a tight 
binding model with nearest and next nearest coupling [25] . 



substrate and N=l properties 


Property 


(jt 


/~i tt 


Uo 


498 K 


195 K 


BA 


24 K 


13 K 


d s 


1.49 A 


1.45 A 


E for 4 He 


-376.15(2) K 


-133.06(1) K 


E for 3 He 


-362.33(1) K 


-126.97(2) K 


BW of 4 He 


9.6 K 


13.6K 


BW of 3 He 


13.7 K 


19.4 K 


m* jra of 4 He 


1.40 


1.05 


m*/m of 3 He 


1.26 


1.01 



Table 5.3: Substrate and N=l properties, with: U ) Depth of potential well; BA) 
Inter-site energy barrier; d s ) Inter-site distance; E ) Ground state energy; BW) 
Bandwidth; m* jm effective mass to bare mass ratio. 



For comparison we have computed with this same method the band energy for 
He on graphite finding substantial agreement with the Carlos and Cole result for the 
lowest band [26J. The bandwidths A of He on these three substrates are given in 
Table IQ1 

From the first energy band it is possible to obtain an estimate of the effective 
mass of one atom of Helium on GF(GH); this is done with a fit of the energy band 



at small wave-vectors (in Fig. |5.9 it is the region around the point T); in fact, for 
small wave-vectors the first energy band ~ h 2 k 2 /(2m*). The effective masses 
m* of the various systems reflect the varying corrugations of the potentials. For 4 He 
( 3 He), the ratios of m* to the bare mass are 1.40 (1.25), 1.10 (1.08) and 1.05 (1.01) on 
GF, graphite and GH, respectively. The smaller mass enhancement of 3 He than 4 He 
reflects the smaller ratio of the corrugation potential to the translational zero-point 
energy. 

5.0.9 Equilibrium density of submonolayer 4 He on GF 

We have studied a 4 He submonolayer on GF. Some of the obtained properties are in 



Tab. 5.4 As He-He interaction we have used an Aziz potential [27]. The ground state 
has been computed for a number of 4 He atoms from 22 to about 120 spanning the 
density range p=0.04-0.09 A~ 2 . On graphite the ground state is the commensurate 
x y/3 R30° state with filling factor 1/3 of the adsorption sites. A similar state 
on GF is obtained by populating fourth neighbor sites (this corresponds to second 
neighbors in one of the sublattices of the honeycomb at a distance 4.482 A) with 
a filling factor of the adsorption sites equal to 1/6 and it corresponds to a density 
P?fe = 0-0574 A -2 . Notice that this density is smaller than the = 0.0636 A~ 2 on 
graphite due to the dilation of the C plane in GF. 
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Figure 5.11: Upper panel: static structure factor for GF at density pfZ. Lower panel: 
the same for GH, at density p™. Lines are guides to the eyes. 



A simple consideration suggests the instability of a similar commensurate state. 
Using the curvature of the He-substrate potential at an adsorption site, the two- 
dimensional zero point energy is estimated to be 55(40) K on GF (GH), much larger 
than the minimum potential barrier 23(13) K, so that such a localized state might be 
unstable. We find indeed that this ordered state is unstable: starting the simulation 
from an ordered configuration after a short Monte Carlo evolution the Bragg peaks 
corresponding to the a/3 x a/3 R30° state disappear and the system evolves into a 
disordered fluid state modulated by the substrate potential. S(k) at this density is 
plotted in Fig. 5.11 as function of k x and k y for two numbers N of particles: the 
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Figure 5.12: Polymer configurations for a liquid system of N = 43 atoms of 4 He 
on GF (left) and iV = 41 atoms of 4 He on GH (right) at equilibrium density Each 
polymer represents the evolution up to 0.4 K -1 in imaginary time of its correspondant 
4 He atom. The substrate is represented in the background with carbon atoms marked 
in gray and the F(H) overlayer marked in red. 

intensity of some of the peaks do not depend on N so they are due to short range 
order, others scale roughly as N and arise from the modulation of the density due to 
the adsorption potential. 



Many— body properties 


Property 


GF 


GH 


Peq 


0.049 A" 2 


0.042 A" 2 


X 


0.142 


0.115 


E 


-377.71(4) K 


-134.02(5) K 


E b 


1.55(6) K 


0.95(6) K 


no 


11 ± 1 % 


22.6 ± 1.3 % 


Ps/p 


0.60(3) 


0.95(3) 


T 

-L C 


0.2-0.3 K 


1.0-1.2 K 


Psat 


0.136 A" 2 


0.108 A- 2 



Table 5.4: Many-body properties, with: p eq ) Equilibrium density; x) coverage; E ) 
Ground state energy per particle; E b ) Binding energy; n ) Condensate fraction; p s /p) 
T=0 K superfluid fraction; T c ) Transition temperature; p sat ) Completion density 



Fig. 5.12 shows a sampled configuration of polymers for a system of 4 He on GF 
and GH at equilibrium density. The spread in space of a single polymer is related to 
the zero point motion, whereas the center of mass of each polymer gives an idea of 
the spatial order of that configuration. As expected, the polymers stay on average 
over the adsorption sites but, apart from the modulation of the external potential, 
there is not spacial order due to the He-He interaction. It is interesting to note that 
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sometimes zero point motion allows a polymer to stretch toward a near adsorption 
site by crossing a saddle point; this is a dynamic in imaginary time that eventually 
leads a polymer to connect with an adjacent one implementing quantum exchanges 
phenomena in a multi connected geometry. On average, such exchanges are more 
frequent in the configuration on the GH substrate rather than that on the GF, this is 
expected because of the strongest confinement given by the GF adsorption potential. 
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Figure 5.13: Panel (a): Energy per particle of 4 He on GF at T = OK. Panel (b): 
Equation of state of 4 He on GH at T = OK. The inset represents a zoom in the 
region around the energy minimum. In both cases, the used particle numbers ranged 
between N = 60 and iV = 120. The dashed line represents a guide to the eye. Circles 
are liquid densities, Squares represent the commensurate 2/7 phase and Triangles are 
incommensurate densities. 



In Fig. 5.0.9 the energy per particle of 4 He on GF and 4 He on GH at the studied 
densities are reported. In the GF case, the energy per particle has a minimum value 
E = -377.71 ± 0.04 K at the density p eq = 0.049 A" 2 . This lies 1.55(6) K below 
the single particle energy, implying that the ground state is a self-bound liquid. In 
the GH case a similar state is obtained, with Eq = —134.02 ± 0.05 K and a binding 
energy per atom of 0.95(6). For comparison, we note that the strictly 2D cohesive 
energy of 4 He [22] is just 0.84 K and the equilibrium density is p = 0.0436 A" 2 . In 
both the cases a liquid phase has been found at least for densities up to filling factors 
x = 1/4 that for the GF case correspond to a density pf^ = 0.0861 A~ 2 , and for the 
GH case to a density pg£ = 0.0912 A" 2 . 

In Fig. 5.0.9 the local density on the x-y plane is shown for 4 He on GF (left) and 
GH (right) at equilibrium density. These local densities clearly reflect the geometry 



88 




(LOS 
O-DB 
O.D4- 
0.02 





Figure 5.14: Local density (in A -2 units) on the x-y plane integrated along the z 
direction of N = 33 atoms of 4 He on GF (left) and N = 41 atoms of 4 He on GH 
(right) at equilibrium density. 



of the adsorption potential shown in Fig. 5.6 The system is energetically allowed to 
stay in a multi-connected space in which adsorption minima are reachable through 
channels that cross a saddle point of the adsorption potential. Note that although 
the geometry in the two cases is the same, the potential barrier above the F (H) 



overlayer is much lower in the GH case (see Fig. 5.7), this produces an higher degree 



of anisotropy in the GF case and reflects in the local density as a non-zero probability 
to occupy an adsorption maximum in the GH case. 
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Figure 5.15: Static structure factor on the x-y plane of iV = 96 atoms of 4 He on GF 
at equilibrium density (upper left), on GH at equilibrium density (upper right), on 
GH at filling factor x = 1/6 (lower left), on GH at x = 1/6 (lower right). 



The static structure factor of 4 He on GF and GH at equilibrium density as well 
as that at filling factor x — 1/6 are shown in Fig. (5.0.9). The sharp peaks reflect the 
density modulation due to the corrugation of the adsorption potential. The crater like 
structure at smaller k represents short range He-He correlations. It can be noticed 
that short range correlations in the GH case are much less anisotropic than that in 
the GF case, this reflects the smaller corrugation of the adsorption potential of GH. 

We have computed the off diagonal one body density matrix pi(r — r'). As can be 
seen in Fig. 5.0.9 p\ reaches a plateau at large r — r' and the Bose Einstein condensate 
(BEC) fraction is 10.3 ±0.4 % at p = 0.049 A~ 2 and 7.3 ±1.5 % at p 1/6 ; the system is 
superfluid. We reach a similar conclusion in the case of the GH substrate: the ground 
state is a liquid with density p eq = 0.045 A -2 and Eq = —134.28 ± 0.02 K per atom 
and the BEC fraction is 22.6 ± 1.3 % at the equilibrium density and 6.8 ± 0.5 % at 
Pi/6 = 0-0608 A -2 . Note that this condensate fraction is significantly smaller than 
the value (~ 40 %) for 4 He in 2D [28]. The smaller value is a consequence of the 
spatial order, albeit imperfect, induced by the substrate potential and of the smaller 
effective surface available to the atoms due to the strong channeling induced by that 
potential. 



In Fig. (5.0.9) the superfluid fraction p 8 J p is shown in function of the temperature 
for a system of iV = 26 atoms of 4 He on GF and N = 20 atoms of 4 He on GH at their 
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Figure 5.16: Panel a) Off diagonal one body density matrix for 4 He on GF at p^ F = 
0.049 A~ 2 (open circles) and p^Z = 0.0574 A -2 (filled circles). Panel b) the same for 
GH, with p™ = 0.042 A -2 and pf£ = 0.0608 A -2 . Lines are guides to the eyes. 



respective equilibrium densities. At finite temperature, the superfluid fraction has 
been estimated with the winding number method. The data at zero temperature has 
been obtained with the evolution of the center of mass of the system for sufficiently 
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Figure 5.17: The superfruid fraction as function of temperature for a system of 
iV = 20 atoms of 4 He on GH and a system of N = 26 atoms of 4 He on GF. The 
transition temperature can be roughly estimated in the range 1.0 - 1.2 K for 4 He- 
GH and 0.2-0.3 K for 4 He-GF. In the inset is displayed the diffusion of the center 
of mass from which the zero-temperature estimation of the superfluidity has been 
extrapolated. 



large imaginary-time [29] r: 



D(t) 



N 



^= lim D{t) 
Rcm (r) - Rcm (0)1 



2d\ 



(5.15) 



where A = h 2 /2m, N is the number of particles, d is the number of dimensions along 
which the contribution to the superfruid fraction is considered, and the squared dis- 
tance [Rchdij) — -Rca/(0)] is evaluated without invoking periodic boundary conditions, 
i.e. including boundary crossing if Rcm{t) leaves the main simulation box. Note that 



in general the estimator for the superfluid fraction in equation |5.15| can be used with 
a PIGS algorithm only when the Hamiltonian of the system explicitly breaks the 
translational symmetry as in the present case. This is a necessary condition because 
even if one starts from a trial state in which the translational symmetry is broken, 
the imaginary-time evolution of \I/t obtained via PIGS (i.e. \1/ T = G(t)^t) imme- 
diately, for every imaginary time r restore the translational invariance unavoidably 
disturbing the estimation of D(t). This is due to th fact that the imaginary time 
evolution depends only on the Hamiltonian via G(t) = exp(-rH). 

It is noticeable that p s j p for He on GH joins smoothly with the T = OK value. 
This is a strong test on our algorithms since these values come from completely 
different computations, however, in the case of GF, the low transition temperature 
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does not allow to reach the T «T C regime. From p s /p at finite T and taking into 
account size effect we estimate the superfluid transition temperature T c ~ 0.2 — 0.3 
K for GF and 1.0 — 1.2 K for GH. From the zero temperature computation of the 
superfluid density we predict that a submonolayer of 4 He film on GF (GH) is an 
anisotropic superfluid with superfluid fraction p s /p = 0.95(3) for GH and 0.60(3) for 
GF. Remarkably, this quantity is less than unity and this is in agreement with the 
predictions by Leggett for a nonuniform superfluid |30j . 
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Figure 5.18: Left: Excitation spectrum, E(k), of 4 He on GF at equilibrium density 
along x and y directions extracted from the position of the quasi-particle excitation 
peaks in the dynamical structure factors obtained via the GIFT algorithm. The 
error-bars represent the 1/2-height widths. The Bijl-Feynman spectrum, Ep(k), is 
also shown. Right: The same for 4 He on GH. Lines are guides to the eye. 
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Dynamics at equilibrium density 

Information about dynamical properties can, in principle, be extracted from imagi- 
nary time correlation functions, without relying on any approximation, focusing on 
an ill-posed inverse problem, i.e. the inversion of the Laplace transform which con- 
nects a suitable imaginary time correlations function /(r) to the relevant spectral 
function. If one consider the dynamical structure factor S(k,u), which is measurable 
in an inelastic neutron scattering experiment, the related imaginary-time correlation 
function is the so called "intermediate scattering function" F(k,r): 

F(k,r) = ^(e TA P % ^ T& P-k) = J duje^S^uj) . (5.16) 

The expression of F(k, Tj) = (e nH p^e~ nH p__j:) I 'N can be estimated via "exact" Quan- 
tum Monte Carlo methods for a discrete set of imaginary time instants Tj. However, 
the extraction of S(k, u) from the above integral equation, based on the limited and 
noisy knowledge of F(k,r), is an ill-posed inverse problem; in fact, the kernel e" WT 
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Figure 5.19: Left: Static structure factors S(k) and strength of the quasi particle peak 
Z(k) as function of q of 4 He on GF at equilibrium density along x and y directions. 
Right: The same for 4 He on GH. 



is strongly smoothing and infinite dynamical structure factors turn out to be com- 
patible with the information on the correlation function, i.e. with F(k, Tj) for the 
different r». Recently, we have developed a technique to face such problems quite 
in general: the Genetic Inversion via Falsification of Theories (GIFT) method[31J. 
GIFT extracts information on spectral functions by averaging among models found 
compatible with observations (i.e. the correlation function for a discrete set of imag- 
inary time instants: /(t»)) via a genetic-algorithm-exploration of a given wide space 
of model spectral functions. When applied to bulk liquid 4 He at T = K, GIFT has 
been found able to extract more information about S(k, u>), separating quantitatively 
the elementary excitation peak from the multiphonon contributions [3"Tj 132] . 

Here we have applied the GIFT algorithm to extract information on excited state 
properties of the equilibrium superfluid phases of 4 He on GF and GH. Intermedi- 
ate scattering functions have been computed for different wave vectors with the 
PIGS method. Well-defined single excitation peaks and multiphonon contribution 
are present in the reconstructed dynamical structure factors via the GIFT algorithm. 
In Fig. 5.18 the position of these peaks as a function of the wave vectors in two 
different directions are shown. In the left panel, which correspond to the GF case, 
the spectrum is found to be highly anisotropic with roton excitations lower than 2 K 
along the x direction and of about 3.5 K along y. This is again a consequence of the 
strong and anisotropic corrugation of the GF substrate respect to the GH case where 
we found a much more isotropic spectrum, with a shallow roton minimum near 5 K. 

Given the extracted single quasi-particle energies of the excitation spectrum one 
can estimate the Landau critical velocity, v c = min(E(k)/hk), for both cases and 
directions; in the GF case, these turn out to be v c ~ 13 m/s along x and v c ~ 31 m/s 
along y, in the GH case we obtain v c ~ 45 m/s along x and v c ~ 51 m/s along y. 
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Figure 5.20: Left: Static density response function, x(™); as function of q of He on 
GF at equilibrium density along x and y directions. Right: The same for 4 He on GH. 



By integrating S(k, u) with respect to u in the range of the sharp peak and in the 
remaining frequency range we can determine the strength of the single quasi-particle 
peak, Z(k), and to the contribution to the static structure factor ? S(k ), coming from 

for both cases, 



5.19 



multiphonon excitations. The results for Z{k) are shown in Fig. 
together with the static structure factor along the same direction; from the ratio 
between Z(k) and S(k) one can measure the efficiency of the single quasi-particle 
excitation channel. 

The efficiency is specially high along the x direction of the GF case where we 
found the roton with the lower energy. 
Also, through the relation 



X (k) = -2p I du 

10 



S(k, uj) 



(5.17) 



one can c ompute the zero temperature static density response function x{k)- In 
Fig. 5.20 we present our results for XW °f 4 H e on GF and GH computed along 



different directions respect to the substrates. 



5.0.10 Properties at high coverages 

The properties of the first layer at high density have been studied. 

At x = 2/7 (p = 0.0984 A" 2 on GF,p = 0.105 A" 2 on GH) on both substrates we 
find that a commensurate triangular solid is stable, or at least metastable, containing 
4 atoms in the unit cell of the triangular lattice rotated by 19.1° with respect to the 
substrate potential. In the unit cell one of the 4 He atoms is localized on an adsorption 
site in the middle of a graphene hexagonal ring, other two atoms approach adsorption 
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Figure 5.21: (Upper panel) Local density (in A -2 units) on the x-y plane of the 2/7 
phase of 4 He on GF compared with the geometry of the substrate. Red balls are 
centered on the position of Fluorine atoms and the green ones on the Carbon atoms. 
Thin white lines enclose the unit cell of the commensurate 2/7 phase. (Lower panel) 
Static structure factor on the k^-k^ plane of the 2/7 phase of iV=112 atoms of 4 He on 
GF. k x and k y axis are expressed in A -1 . Red arrows point to the peaks corresponding 
to the density modulation imposed by the adsorption potential. 



sites of the other kind and finally the fourth one is centered on a saddle point of the 
potential. This state has some similarity with the 4/7 commensurate state found 
for 3 He in the second layer on graphite[33]. The local density (Fig. 5.21 for 4 He on 
GF) displays the presence of a superlattice with four atoms in the unit cell of the 
triangular lattice. The static structure factor S(k x , k y ) has prominent Bragg peaks 
forming three stars. S(k x , k y ) for 4 He on GF is shown in Fig. 5.21 The star of the six 
highest peaks is the one of a triangular lattice with lattice parameter equal to that of 



96 



>'!*> 




Figure 5.22: (Upper panel) Local density (in A -2 units) on the x-y plane of the 
2/7 phase of 4 He on GH compared with the geometry of the substrate. Red balls 
are centered on the position of Hydrogen atoms and the green ones on the Carbon 
atoms. Thin white lines enclose the unit cell of the commensurate 2/7 phase. (Lower 
panel) Static structure factor on the k x -k y plane of the 2/7 phase of iV=112 atoms 
of 4 He on GH. k x and k y axis are expressed in A -1 . Red arrows point to the peaks 
corresponding to the density modulation imposed by the adsorption potential. 



a triangular lattice at this density. Another star represents the density modulation 
due to the adsorption potential. The third star formed by six less intense peaks 
at a smaller wave vector is a combination of vectors of the two previous stars, thus 
corresponding to interference between the triangular and the honeycomb modulation. 
The intensity of all these peaks scale with the number of particles (data not shown). 
Additional peaks are present reflecting the superlattice but these peaks are very weak 
and hardly visible in the figure. Returning to the local density in Fig. 5.21 it can 
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Figure 5.23: Polymer configurations for a system of iV = 56 atoms of 4 He on GF (left) 
and GH (right) at filling x = 2/7. Each polymer represents the evolution up to 0.4 
Kr 1 in imaginary time of its correspondent 4 He atom. The substrate is represented 
in the background with carbon atoms marked in gray and the F(H) overlayer marked 
in red. 



be noticed that some of the spots, those located at a saddle point, are elongated 
indicating that the atoms visit also the neighboring adsorption sites. S(k x , k y ) and 
the local density of 4 He at coverage 2/7 on GH are shown in Fig. 5.22 The results 
are similar to those on GF, it can be noticed the much smaller intensity of the Bragg 
peaks due to the adsorption potential in the case of GH as it can be expected due to 
the weaker corrugation of the adsorption potential of GH. 

It should be noted that in this 2/7 state not all atoms are localized around a single 
adsorption site but some atoms visit two or three neighboring sites, as consequence 
there is spatial order but the atoms are rather mobile and exchange easily so these 
solids might be supersolid. Evidence of this is indeed what we find at T = OK for both 
the substrates with the superfluidity estimation through the diffusion of the center 
of mass of the system in imaginary time (see Fig. 5.25). At the commensurate 2/7 
phase we estimate a superfluid fraction of 0.23 for GF and of 0.61 for GH. 

At coverages around 2/7 we find that 4 He has an incommensurate triangular 
order deformed by the substrate potential and defected because such order is not 
compatible with the periodic boundary contitions at the box sides. We discuss first 
4 He on GF. We have investigated the density range between p^Z = 0.0984 A -2 and 
pfat = 0.136 A~ 2 and as an example S(k x , k y ) at p = 0.123 A -2 is shown in Fig. 5.27 
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Figure 5.24: Polymer configurations in the incommensurate solid density range for a 
system of N = 82 atoms of 4 He on GF at density p = 0.118 A~2 (left) and a system 
of TV = 98 atoms of 4 He on GH at density p = 0.0916 A _ 2 (right). Each polymer 
represents the evolution up to 0.4 K -1 in imaginary time of its correspondent 4 He 
atom. The substrate is represented in the background with carbon atoms marked in 
gray and the F(H) overlayer marked in red. 



and 5.26 As initial configuration we have used a disordered one as well as an or- 
dered triangular configuration. The runs converge to the same results: S(k x ,k y ) is 
dominated by a star of six peaks as expected for a triangular solid. The wave vectors 
of these peaks are not exactly equal to the value of an ideal triangular solid at this 
density implying that the triangular order is deformed in order to better fit within 
the simulation box. S(k x , k y ) has additional Bragg peaks corresponding to the mod- 
ulation of the substrate potential and to the interference between the previous two 
sets of peaks. Additional smaller peaks are present presumably as consequence of the 
presence of defects. The modulus of the main Bragg peaks increases in a smooth way 
as the density is increased as expected for a triangular solid. The observed deviations 

1/2 

from the value ks = 47r (p/2y / 3) of an ideal triangular solid are explained by the 
deformations of the lattice and by the presence of some defects, mainly dislocations, 
that can be observed from the configuration of the atoms (data not shown). 
In the case of GH we have investigated the density range 0.0916 A~ 2 - 

0.105 A~ 2 . Again S(k x ,k y ) is dominated by the Bragg peaks of a triangular lat- 
for S(k) at density p 



5.27 



GH 

P2/7 



0.102 A ) that is incommensurate with 



tice (see Fig. 
respect to the substrate periodicity. 

In Fig. (5.23 ) and (5.24 we show a sampled configuration of polymers for a system 
of 4 He on GF and GH at respectively commensurate and incommensurate density. 
In the commensurate density case, Fig. 5.23 the adsorption potential has a greater 
influence in the GF case and causes the polymers on the saddle point to spread to 
neighboring adsorption sites; this effect is much less evident in the GH case where 
the shape of the polymer is more isotropic and the occupation of adsorption maxima 
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Figure 5.25: Center of mass diffusion for the 2/7 phase of GF (Circles) and GH 
(Squares) and the 1/3 phase of Graphite (Triangles). The superfTuid fractions are 
obtained as the long-r limit of the plotted functions. The horizontal dashed lines 
represent the value of the superfluid fraction, 0.23 for GF and 0.61 for GH. The dashed 
exponential curve is a fit to Graphite data and has a vanishing long r behavior. 



on top of the H overlayer is more likely. In the incommensurate case, Fig. 5.24, the 
presence of point dislocations are clearly visible in both the cases. 

The static structure factors in the solid phase show a characteristic structure of 



three sets of six peaks that is represented in figure (5.26). The set that has the 



higher intensity represents the Bragg peaks of a triangular lattice. The six peaks 



that in figures 5.21 , 5.22 and 5.27 are marked by the red arrows represent the density 
modulation induced by the adsorption potential like the peaks in Fig. 5.0.9 The 
third set of peaks is merely an interference pattern of the first two sets. 

Increasing the number N, at some point some atoms spill out of the first layer and 
the density profile in the direction normal to the surface develops two well separated 
peaks. We have thus estimated the first layer's completion density, pf^ ■ The 
promotion to the second layer takes place at a density pf£ = 0.136 A~ 2 for the 
GF case and a density p^t — 0.108 A~ 2 for the GH case. Beyond such densities, 
the occupation of the second layer is clearly visible as a secondary peak in the local 



density along the ^-direction displayed in Fig. 5.0.10 



5.0.11 Equation of state of 3 He on GF and GH 

The ground state of 3 He on graphite is the v^3 x y/3 R30° state. We expect that the 
analogous commensurate state on GF and GH is unstable, as for 4 He, because the 
smaller mass makes 3 He localization more expensive. 
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Figure 5.26: Level curves of the static structure factor of the iV=86 incommensurate 
solid phase of 4 He on GF (Figure 5.27). The six peaks marked (A) reflect the density 



modulation due to the corrugation of the adsorption potential. The six peaks marked 
(B) represent the periodicity of a triangular lattice and the peaks marked (C) are 
interference patterns from the density modulations represented by (A) and (B). k x 
and k y axis are expressed in A -1 . 



The ground state energy as function of density of mass 3 bosons and of the 
Fermionic 3 He on GF and on GH are plotted in Fig. 5.29 as function of density. In 



both cases, the system under study was composed of N =18 atoms of 3 He. 
In the figure we plot also the 3 He energy based on the crude approximation of taking 
as Fermi-Bose gap the kinetic energy Kf ree = h 2 np/2m* of free fermions, where m* is 
the effective mass of a 3 He atom on the substrate (m*/m = 1.26 for GF, m*/m = 1.01 
for GH). 



As shown in Fig. ( J5.29 ) , the a/3 x a/3 R30° commensurate state for a mass 3 
boson system is indeed unstable toward a fluid state on both substrates, in fact, the 
energy at the density corresponding to the v3 x a/3 state is well above the energy 
at lower densities implying that this ordered commensurate state is indeed unstable 
and the system is in a fluid state. As a consequence we predict the existence of two 
new anisotropic Fermi fluids, in the sense that the local density is non-uniform and 
anisotropic, with a tunable density depending on the 3 He coverage. The density range 
depends on whether the 3 He atoms form a self-bound state. Such a self-bound state 
seems unlikely to occur for 3 He on GH on the basis of our computations. On the 
contrary a self-bound state might be present on GF. For mass 3 Bosons we find a 
bound state with a binding energy E Q = —0.22 K at density p eq = 0.03 A" 2 . Adding 
to the boson energy the Fermi-Bose gap the energy yields a shallow minimum in the 
density range 0.015-0.025 A" 2 . The energy per particle at this minimum is equal 
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Figure 5.27: (Left) Static structure factor on the k x -k y plane of the incommensurate 
solid phase of iV=86 atoms of 4 He on GF at density p =0.123 A -2 . k x and k y axis 
are expressed in A -1 . Red arrows point to the peaks corresponding to the density 
modulation imposed by the adsorption potential. (Right) Static structure factor on 
the k x -k y plane of the incommensurate solid phase of iV=66 atoms of 4 He on GH at 
density p =0.102 A~ 2 . k x and axis are expressed in A -1 . Red arrows point to the 
peaks corresponding to the density modulation imposed by the adsorption potential. 



within the statistical error to the energy of a single adsorbed 3 He on GF so that the 
existence of a self-bound state on GF is an unresolved possibility. 

Remark: An accurate approximation for the energy per particle for 3 He on GF 
and GH has been obtained via the Fermionic Correlations technique [Ml [35] . This 
methodology has been explained in Chapter [3] given a specific Hamiltonian, the 
Fermionic Correlations technique extracts the energy gap between the symmetric 
and antisymmetric ground state from a suitable Fermionic imaginary-time correlation 
function computed as an exact average on the Bose ground state: 

Cf{t) = V TJBTTW- > ( 5 ' 18 ) 

where Af is, typically, a Slater determinant. The lowest energy contribution in 
Cf{t) yields the exact gap between the Fermionic and the Bosonic ground states, 
provided that one is able to obtain the inverse Laplace transform of Cf{t); this can 



be readily seen by formally expressing (5.18) on the basis {ip^ } n >o of eigenvectors of 



H corresponding to the eigenvalues {E^} n 



>o- 



• I(-4f«A b IV'„ f >« ( a.)I 



We have shown that this analytic continuation procedure can be handled effi- 
ciently with statistical inversion procedures, like the GIFT algorithm introduced in 
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Figure 5.28: Local density along the z-direction of 4 He on GF (with N=lll) and 
4 He on GH (with N=79) at a density beyond the promotion density. The occupation 
of the first and the second layer are clearly visible as two peaks. The area under the 
peaks represents the number of 4 He atoms in the corresponding layer. 



Ref.[3T].The Fermi-Bose gap Eq - is an extensive quantity, so this method can 
be applied provided that the system is not too large. 

5.0.12 Discussion 

He adsorption on new substrate materials is valuable because of the fundamental 
importance of helium in many-body physics, with a variety of phases seen in both 
2D and 3D. Our results indicate that the GF substrate provides the strongest binding 
of any surface (since the previous record was held by graphite). Moreover, the novel 
symmetry, the smaller intersite distance and large corrugation imply that quite novel 
properties may be anticipated for this system. This is indeed the case. When many 
4 He atoms are adsorbed on GF and on GH a very striking result is that the ground 
state is a low density liquid modulated by the substrate potential and the system has 
BEC, i.e. it is a superfluid. This is qualitatively different from graphite for which 
the lowest energy state is the \/3 x R30° commensurate one with no BEC [36J. 
We have verified that such an ordered state on GF and GH is unstable relative to 
the liquid phase. It should be noticed that some of the parameters in the adsorption 
potential are not known with high precision or they have been adopted from other 
systems. We have verified that even a change of parameters like a, Cqf and /3 by 20% 
does not modify the qualitative behavior of the adsorbed He atoms even if there can 
be a sizable change in the value, for instance, of the adsorption energy. Measurement 
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Figure 5.29: Panel a). Ground state energy as function of density of mass 3 bosons 
on GF (circles), Fermionic 3 He-GF obtained via the Fermionic correlations method 
(triangles) and Fermionic 3 He-GF obtained by approximating the Bose-Fermi gap 
with the kinetic energy of the free fermion gas (diamonds). Panel b). Same as for 
panel (a) for 3 He on GH. 



of thermodynamic properties and He atomic beam scattering experiments from GF 
and GH will be important to test the accuracy of our model potentials. A remarkable 
result is the superfluid behavior of the 2/7 phase that, however, might be a property 
of the system at strictly T = OK and is non reachable by experiments; on the other 
hand there might truly be a "supersolid" phase transition at a temperature in the 
mK range that is not accessible by QMC computations. This, together with all the 
novel phenomena for He atoms on GF and GH that have been predicted in this 
work, calls for experimental verification. There is also an important aspect that 
should be considered in view of experiments; it might be difficult to have a 100% 
reacted graphene sheet with fluorine. However, the presence on GF of small regions 
of unreacted graphane should not affect the properties of the adsorbed film because 
the He or H2 atoms are preferentially adsorbed on the F covered regions of graphene. 
This behavior, however, may change when coverages beyond the first layer completion 
on GF and GH are considered; in such cases it could be that the adsorbed atoms begin 
to populate the unreacted regions. The QMC techniques here employed may be used 
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to investigate also these interesting cases given that a suitable interaction potential 
is provided. 

From the theoretical point of view many extensions of the present computations 
can be foreseen, for instance the characterization of the commensurate 2/7 phases on 
GF and GH, of the system under rotation and the study the phase diagram of p-H 2 
on GF. As a perspective of future work we plan to provide predictions concerning the 
phase diagrams and thermodynamic properties for both He/GF and He/GH, hoping 
to stimulate experimental studies of these systems. 
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Chapter 6 
Conclusions 



The idea underlying this work has been the study of strongly interacting quantum sys- 
tems along with the development of new methodologies in the field of QMC. Strongly 
interacting quantum systems are indeed a fascinating field of research, much is yet 
unknown and indeed a proof of this assertion might be found in our results: we have 
studied new adsorbed phases of 4 He and also predicted the presence of a modulated 
superfluid given by the interplay between interatomic potentials and quantum tun- 
neling. Strongly interacting Fermi systems are even more unexplored due to the sign 
problem. The methodological aspect of this work has been thus focused to develop 
a technique that can study the dynamics of such systems. This technique is an evo- 
lution of the Fermionic Correlations method; we have shown that, even though this 
methodology becomes unpractical for big numbers of particles, it indeed can compute 
an ab-initio low-energy excitation spectrum of two-dimensional 3 He. 

In the conclusions of this work, we remind the main results obtained and presented 
throughout this work; we have already drawn conclusions in each chapter, here we 
will comment mostly the computations that are still in progress and even the "failed 
attempts" . 

2d 3 He. Our simulation of two-dimensional 3 He gave a spin susceptibility as func- 
tion of density that is in very good agreement with experimental data; our obtained 
polarization curves indicates that the ferromagnetic fluid is never stable and the sys- 
tem crystallizes into a triangular lattice from the paramagnetic fluid at a density 
of 0.061 A -2 . With an extension of the Fermionic Correlation (FC) technique, we 
have been able to obtain the first ab-initio evaluation of the zero-sound mode and 
the dynamic structure factor of 2c? 3 He that is in remarkably good agreement with 
experiments. This excitation spectrum, moreover, turned out to have striking simi- 
larities with the phonon-maxon-roton spectrum of 4 He; this indicates that the effects 
of the inter-atomic potential, in particular its strong repulsive part, dominate over 
the effects of the quantum symmetry. Another interesting question is whether the 
zero-sound mode, which is known to enter the particle-hole band, reemerges at wave 
vectors corresponding to the "roton" minimum of the spectrum; it is possible, with 
the FC method, to compute the particle-hole excitations and indeed we have pre- 
sented preliminary results that show that the "roton" is still in the particle-hole 
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band. However, the re-emergence of the roton is a rather difficult question to answer 
with a simulation of a finite system. This is because there are relevant size effects 
on the particle-hole; these effects are due to the fact that we are far from the ther- 
modynamic limit and the particle-hole is not a continuum; in order to obtain more 
conclusive data, a scaling analysis on bigger systems is in order. Such a study that 
is very demanding in term of computing resources and has been planned for future 
work. 

We also attempted a study of the spin-waves excitations but in this case we found 
that the results were highly dependent on the direction of the wave vector. This 
anisotropy is a clear sign of size effects; unfortunately, these size effects for spin- 
waves are stronger than in the zero-sound or even the particle hole case and thus 
require the study of systems with particle numbers for which, like for most QMC 
methods, FC becomes unpractical. 

An interesting perspective is the application of the FC technique to the study of 
elementary excitations of the 2d electron gas; possibly this would provide an ab-initio 
evaluation of the plasmon excitations. 

4 He on Graphene-Fluoryde and Graphane. The study of 4 He adsorption on 
Graphene-Fluoryde (GF) and Graphane (GH) has been a comprehensive and artic- 
ulated work. At the early stages of the project we showed that the commensurate 
x R30° phase is unstable on both substrates. We then determined the equilib- 
rium density at T = K; our results indicated clearly that on both the substrates the 
equilibrium density has a condensate fraction and is thus a modulated superfluid. We 
determined the superfluid fraction at zero and finite temperature giving also a rough 
estimate of the fluid-superfluid transition temperature. The study of the equilibrium 
density at T = K comprised also the excitation spectrum, and we have shown the 
phono-roton spectrum of 4 He on GF and GH. We focused then on high coverages of 
the monolayer and found a density range, not yet precisely determined, in which 4 He 
forms an incommensurate triangular solid; a remarkable result is that on both GF and 
GH a commensurate phase at filling factor x = 2/7 is stable or at least metastable. 
This result becomes even more interesting because we found a first evidence of super- 
fluidity at zero temperature: at this density, the system may possibly be both solid 
and superfluid, in other words this system could posses the long sought property of 
supersolidity. 

For the immediate future, we plan to study further this commensurate density with 
also a size scaling aimed to better estimate the finite size effects on the properties of 
the 2/7 phase, in particular on the superfluid fraction. In our further studies, there 
are mainly two points to inspect: first, is the 2/7 phase thermodynamically stable? 
This far, we have shown that it is mechanically stable, meaning that, at the density 
corresponding to x — 2/7, the configuration that gives the lowest energy is the 2/7 
triangular lattice; now we are planning to search for signatures of possible phase 
transitions between the incommensurate solid and the 2/7 phase. Second, we are 
searching for more evidence of superfluidity in the commensurate phase, we already 
tried the computation of the superfluid density at finite temperature down to 0.5 K but 
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we did not find any superfluid signal, indeed, the very low rate of exchanges between 
atoms suggests that, if any, the transition temperature to the supersolid state would 
be too low to be reachable with PIMC with our current computing resources; another 
approach is the computation of the one body density matrix at zero temperature, in 
fact, the presence of even a small fraction of condensate would be a strong support for 
the supersolidity of this phase; this is very demanding in term of computing resources 
and is planned for the next future. 

Besides the supersolidity, we have scheduled also a deeper study of the incom- 
mensurate density range with a quantitative characterization of the defects. This will 
prepare the background for the study of the second adsorbed layer that will be left 
for future works. 

We conclude with a last remark. For the Helium-substrate interaction we have 
adopted a semi-empirical potential: its repulsive part has been obtained from a DFT 
calculation of the electron density of the substrate, on the other hand, the attractive 
part has been modeled with a Van der Waals type interaction with parameters taken 
from literature, adopted from the interaction potentials of Helium with similar chem- 
ical compounds. This study, as consequence, can be considered a semi-quantitative 
approach but the very fact that we find qualitatively the same behavior on two sub- 
strates that have completely different values of energies and corrugation is a strong 
proof of plausibility of our results; moreover, the robustness of the results has been 
explicitly tested at the y/3 x \/3 R30° density with a variation up to 20% of the 
parameters of the He-substrate potential. More exactly, this is an "exact" study on 
a semi-empirical Hamiltonian aimed to the research of new properties of adsorbed 
matter, our hope is that this predictive work will encourage, on one side, the develop- 
ment of more accurate Helium-substrate interaction potential and, on the other side, 
the experimental exploration of this subject so fascinating and full of surprises. 
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Chapter 7 

Computational details 



In this Chapter some technical details of the Monte Carlo techniques introduced 
in Sec. [2] will be described. Monte Carlo sampling will then be applied to the 
problem of evaluating physical properties of quantum systems at both zero and finite 
temperature. 

The basic idea underlying the used path integral methods is that the computation 
of an expectation value in a quantum system can be viewed as an iV-dimensional 
integral; in the case of a bosonic system, the ground state wave function can be 
chosen real and non-negative and this integral can be interpreted as the average of a 
random variable over a probability density [3J. 



7.1 Monte Carlo integration: the strategy 

An effective way to compute an iV-dimensional integral is to employ Monte Carlo 
(MC). Monte Carlo basically means "the use of random numbers in order to solve 
a problem". In our case, the problem is the iV-dimensional integral representing an 
expectation value for a quantum many-body system. In order to show how MC is 
employed in our context, consider, as an illustrative example, a function f(x) that is 
a product of an arbitrary function g(x) and a probability density p(x) 

f(x) = g(x) ■ p(x) with (7.1) 
p(x) > OVf G r, / dxp(x) = 1 (7.2) 



The integral can be rewritten as[l], 

1 N 

/ dxg{x)p{x) = Jim ^7 Y%0?i) ^ (g)j 

JL i=l 



(7.3) 



where Xi are elements sampled from the probability density p(x). The integral value 
is thus the average of g(x) over sets of values x that are sampled from the probability 
density p(x). The advantage of MC is that, given a way to sample p(x), the computing 
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time required for the evaluation of the integral does not scale with the dimensionality; 
this is very important since, as will be clear further on, the integrals that are computed 
in this context have generally a very high dimensionality. 




Figure 7.1: a) Schematic ID representation of the sampling of a probability density 
and its discretization in N bins, b) The partition of the unity used in the faked 
roulette method: the interval [0; 1) is divided into N bins, the m bin is the interval 



Yl^JLi P( x j)/^'i Y^JLi p( x i) / J > where the normalization Z = Yl!j=iP( x j)- The 
faked roulette is a method to randomly chose m: a random number r uniformly 
distributed in [0; 1) is generated, m corresponds to the interval I m in which r falls, 
namely: I m so that I m fl {r} = {r} 



This however requires that one is able to sample an arbitrary iV-dimensional prob- 
ability distribution. Sampling means the generation of a random variable according 
to an arbitrary probability density p(x); as sketched in Fig. 7.1 , sampling can be done 
by dividing the domain of p{x) in K N bins with an assigned probability p\ = 
with Z a normalization constant and 5?,- the central coordinate of the z-th bin. A sim- 



ple way to extract a random variable value is through a faked roulette (see Fig. 7.1), 
however, the computational weight of this approach increases exponentially with the 
number of degrees of freedom and is thus unpractical for the evaluation of Eq. (7.1). 
A more sophisticated and efficient way to sample an arbitrary probability density is 
with Markov chains[5j. As is shown in the next section, a Markov chain has at least 
one invariant probability density and there is a sufficient condition for its unique- 
ness; the basic idea is thus to build a Markov chain that converges to the required 
unique invariant probability density. In the next section we provide a mathematical 
demonstration of the properties of the Markov chains used in this context, after that 
section an algorithm that can be used to build the required Markov chain, namely 
the Metropolis algorithm |6J, will be described. 
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7.1.1 Mathematics of Markov chains 



In this section, we follow Ref. [2] and show the mathematical basis of the Markov 
chains. Let's consider a given finite set E = {1, ...N} and a probability space 
(O, J 7 , P), where Q is a sample space, J 7 is a a-algebra on Q and P a probability 
measure. 

Definition 1. A Markov chain on a sample space E is a sequence of random 
variables {X n } neN o, X n : VL ^ E such that there are non negative numbers Vi^.j(n), 
n G N° and i,j&E for which, 

P( X n+i = j\X n = h X n _i = i n -\, ...,X = i ) = 

P{X n+1 = j\X n =i) = V^ 3 {n) (7.4) 

whenever the conditional probabilities P(-\-, ■■■) are defined. From here on we focus 
on Markov chain that are independent on time translations, in this case (n) does 
not depend on n. 

Transition matrix. The non-negative numbers Vi-^j can be represented in an 
N x N matrix V_ that is referred as transition matrix of the Markov chain. Following 
from the definition of conditional probability, the sum of the element of a row is one, 
namely = 1 for i = 1,...N. The probability distribution of the random 

variable X is the starting probability of the chain and is defined by the numbers 
Vk = P(X = k) , k G E, this probability distribution can be identified by a row 
vector in H N , v = (vi, ...,Vn). 

Statement: a Markov chain that is independent on time translations is defined 
by a starting probability and a transition matrix. This can be shown if, from the 
definition of conditional probability, we first obtain the probability distribution of 
X x : 

N N 

P{X X = k) = ^P(X = h)P(X! = k\X = h) = Y,VhVh^ k (7.5) 

h=l h=l 

that in matrix notation becomes, 

v {1) = vV (7.6) 

where we define vj 1 ^ = P{X\ = k). Iterating to obtain the probability distribution 
for the next time step, P(X 2 ), 

N 

P{X 2 = k) = J2 p ( x i = l ) p ( x 2 = k\X! = l) = 
1=1 

N N N 

= p ( x i = i)p^k = J2J2 ^-+1^ = 

1=1 1=1 h=l 

N N 

h=l 1=1 
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in matrix form: 



= ^p2 ^g) 

iterating this rule the probability distribution at time step n can be obtained as the 
n-power of the transition matrix applied to the starting probability, = vV_ n . 
It will be used later on in this section the m steps transition probability, Vj;™] = 
P(X n+m = j\X n = i). 

Statement: v£] are the matrix elements ofV m . This, again, is shown by 
iteration: 

p/y _ -I V _ -\ _ P\X n +m = j, X n = i) 

P(X n = i) 
_ v - ^ P(X n+m = j, X n+m _i = h, X n — i) 

ft T 

P(X n+m = j, X n+m _i = h, X n = i) P(X n+m -i = h, X n = i) 
P(X n+m - 1 = h,X n = i) P(X n = i) 

P(X n+m = j|X n+m _i = h,X n = i)P{X n+rn ^i = h\X n = i) = 



E 

h=i 

N 



h=l 

N 



J2'Ph^ j P(X n+m _ 1 = h\X n = i). (7.9) 



h=l 

Repeating this passage with P(X n+m _ 1 = h\X n = i), and so on until P(X n = k\X n = 
i), the statement is demonstrated. 

Invariant probabilities. Given a probability distribution re on the set E identi- 
fied by a row vector n = (tti, 7Tjv) £ R w and a Markov chain defined by a transition 
matrix V and starting probability v, it is invariant if: 

7L = 7LV, (7.10) 

in the particular case of v = 7r we say that the Markov chain is stationary. 

Theorem 1. (Markov-Kakutani) There is always at least one invariant prob- 
ability distribution. 

Proof. The probabilities on E are mapped onto the set 

S = jx G : < Xi < 1 , J2 x i = 1 J > ( 7 - n ) 

this is a closed and limited set in and hence, by the Bolzano- Weierstrass theorem, 
it is a compact set; given a sequence of points in S, it is thus possible to define a 
subsequence that converges to a point in S. From a point x of S, we define the 
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sequence: 

n-l 



n 

k=0 



this vector has non negative components because is a product of elements with non 
negative components; x n also belongs to S, this is readily seen with a summation of 
its components x n ^ 



N „ n-l AT AT „ n-l AT 



E = \ E E E = ^ E E = 1 ( 7 - 13 ) 

i=l k=0 h=l i=l fc=0 7i=l 

where in the last passage the property of the transition probability £V 7-^™$ = 1 was 
used. Having proved that x n G 5, there is a subsequence {x nfe } converging to a point 
7T G S. Now we write 





-n k -n k 2 — 


n^ — 1 






- E = 


. h=0 


7i=0 / 



= —(x-xV nk ) (7.14) 

taking the limit for k to infinity, we observe that while n k diverges, the quantity 
x — xV nk remains finite because is the difference between two elements of the limited 
set S. We thus obtain the following relation for 7r: 

7l - 7l V= lim (x„ -x„ V)= lim — (x - xV nk ) = 0, (7.15) 

that demonstrates the theorem. □ 

The invariant probability distributions of a Markov chain can be obtained from 
the solution of the linear system 

N 

TTj = ^ITiVi^j, (7.16) 
i=l 

however, a sufficient condition for n to be invariant is that it satisfies the detailed 
balance equation: 

-i'P; .j ;~ iP, ., (7.17) 
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for every i,j G E. The demonstration of this statement comes directly from the 
definition of transition matrix 

n n 

KiVi^j = ^ KjVj^i = TTj (7. 18) 

i=\ i=l 

The detailed balance condition will be used in the next section when the Metropolis 
algorithm will be described. This algorithm is used to build a Markov chain that 
converges to an arbitrary invariant probability density. 

Uniqueness of invariant probability distributions. Markov chains wouldn't 
be so much useful in Monte Carlo if there were not conditions of uniqueness of the 
invariant probability distribution. The uniqueness property is in fact what guarantees 
that the Metropolis algorithm converges to the wanted probability distribution. To 
state and prove this property some definitions are necessary. 

Definition 2. Given the transition matrix V = {Vi-^j}i t j of a time invariant 
Markov chain, 

• V is irreducible if for each i, j G E there is a positive integer number m — m(i, j) 
so that V^] > 0. 

• V is regular if there is a positive integer number m for which > for every 
i,j G E. 

Clearly, a regular transition matrix is also irreducible but the opposite is not generally 
true, however if an irreducible transition matrix satisfies the following criterion, then 
we will show that it is also regular. 

Statement. If a transition matrix is irreducible and there is h G E such that 
'Ph^h > 0, then that transition matrix is also regular. 

Proof. From the definition of irreducibility, for each i,j G E there is m = m(i,j) > 
such that Vj;™] > 0; defined s = maxjj eE m(i,j), then Vji^l > for each I, k G E. 
In fact, one can always use iteratively the transition element Vh^h > to express 
Vi^f. as a chain of products: the irreducibility guarantees that given two elements 
l,k G E there are positive integer n\ = n(l, h) and n2 = n(h, k) such that > 

and Vj^l > 0; the element Vj^l will thus be expressed as: 

<T->( 2s ) ^ ^(n^n^^h))^ ^ ,p(n2=n(/i,£!)) „ (7 -[Q\ 

' l^rk — ' l^rh ' h—th---' h—th' h-^rk > u > l'- iy J 

and this proves the statement. □ 
We state now the uniqueness theorem. 

Theorem 2. (Markov) If a transition matrix is regular, then there is only one 
invariant probability ir and, for any starting probability v, the following holds 

7r J= lim (vV n ), (7.20) 
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Proof. From Markov-Kakutani theorem at least one invariant probability ir exist and, 
by definition 

1L = 1lV, £V fe = l, 0<vr fc <l (7.21) 

k 

Let's consider the one dimension vector subspace of generated by n: 

K = {u G C N \u = tvr, t G C} (7.22) 

By the hypothesis and this definition follows respectively that 7r is an eigenvector of 
V_ and 7[ G V T . Define also the subspace Vo: 



V = Le C N \J2 
{ k 



y k = o\. (7.23) 



This subspace has dimension M — 1 and Vo PI V n = {0} because V^ is made of 
elements u for which = t. Following from this conditions, the vector space 

decomposes in the direct sum Vo © V n . This implies that any element v E C N can 
be uniquely written as: 

v = tK + y, y e V (7.24) 

The eigenvalues equation for V is: 

vV = \v, A G C (7.25) 

Let's consider also that given an element y G Vo, then also the element yV G Vo, 
this is seen in this passage: 

= £Z>n-* = x>I>^* = I> = ( 7 - 26 ) 

fe k i i k i 

N 



For an eigenvector v G C of V, using Eq. (7.24) 



vV = tnV + yV = \tn + \y (7.27) 

here, yV_ £ Vo, and thus, for the decomposition of the vector space C", it must 
necessarily be 

' tnV = Xtn 

Let's consider the eigenvalue equation 

yV = \y, (7.29) 
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with y G Vo and y 7^ 0. Explicitating the matrix product and focusing on the i th 
component, 



(7.30) 



Taking the absolute value, summing over the components, and considering that V_ 
has non-negative elements. 



wEw=E 



that implies 



IAI < 1 



(7.31) 



(7.32) 



If the transition matrix V_ has the elements strictly positive, then Eq. (7.31 ) is strict 



inequality; this can be understood if we consider a summation of complex numbers 



E 



N 



,e^ n , it is true that 



A? 



N 

E 

i=l 



i<t>n I 



(7.33) 



if and only if n = a G R Vn G [1, N] ; however, by hypothesis E^ fc 7/^ = and y ^ 
and this necessarily implies that at least one component of y must have a different 
phase; moreover, we have assumed that the matrix elements of Vk^i are strictly 
positive so in a product they won't change the phases of the vector, meaning that 
ykPk^i has the same phase of y^. Hence, under the condition of strict positiveness of 
the matrix elements of V, holds that 



|A|<1 



(7.34) 



This is the crucial point of this demonstration. Using now the regularity condition 
we demonstrate the following statement. 

Statement. For a regular transition matrix, |A| < 1. 

By definition of regularity, if V_ is regular, then there is an integer m > such 
that V m has strictly positive elements, V m is also a transition matrix, so the passages 
of the demonstration can be applied also to V m resulting in 



|A m | < 1 



(7.35) 



and this immediately implies that |A| < 1. Consider now an arbitrary initial proba- 
bility v, then 



vV n = (K + v-7L)V n 



7T 



zi)E n 



(7.36) 
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From Eq. (7.24) it is clear that v — jiE Vq, but the linear operator V_ defined on Vo 
has eigenvalues that in modulus are strictly lesser than 1; from Functional Analisys[8], 
given these conditions, it follows that 

lim {{v-7 L )V n } = 0. (7.37) 

In conclusion we have, 

lim vP n = n+ lim {(v - n)V n } = vr, (7.38) 

n— >+oo n— >+oo 

this proves the theorem of uniqueness. □ 
7.1.2 The Metropolis algorithm 

Given a probability distribution it defined on a finite set E = {1, N} we now show 
a recipe to obtain a transition matrix V for which ir is the only invariant probability, 
namely 

7T,- = lim [vV n ) j . (7.39) 

n— >+oo 

The recipe that we are going to show is the Metropolis algorithm[6j and it allows to 
sample any arbitrary probability density that satisfies some conditions. This algo- 
rithm is relevant in our context because it is used to evaluate iV-dimensional integrals 



like the one in Eq. (7.3) 



Theorem 3. (Metropolis). Given a strictly positive probability distribution ty, 
7Tj > V j ; that is not the uniform probability density, for each probability distribution 
v, there is a Markov chain with initial probability v and regular transition matrix V 
that has it as invariant distribution probability. The transition matrix V_ is defined 
as: 

{A v ' ^ ./• TTj > 7Ti 
i±h ^<Ti t (7.40) 

where £ is any symmetric and irreducible transition matrix. 

Proof. With this choice of V_, let's show that it satisfies the detailed balance condition; 



chose two elements i,j of E such that Tij < m, applying equation (7.40) we obtain 



KiPi^j = KiCi-^j— = Ci->jTTj = 7CjVi->j, (7-41) 

where we used the symmetry of £; this shows that it is indeed an invariant probability 
associated to V_. To prove the theorem we have to show that tt is also the unique 
invariant probability of V; to this purpose, we will show then that V is regular. Let's 
start from the irreducibility. 
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Statement. V is irreducible. This follows directly from the irreducibility of C. 
The irreducibility is in fact related to the non-zero element of the transition matrix, 
this implies that if a transition matrix is irreducible, then it is also true for any 
transition matrix that has at least the same non-zero elements. This is exactly the 



case for V as can be readily seen from its definition in Eq. (7.40). To show that V_ is 
also regular, as shown in the previous section, it is enough to verify that there is an 
element i G E such that Pj ^i > 0. By hypothesis, tc is not the uniform probability 
distribution and thus there is a subset M C E in which tt is maximum, moreover, 
due to the irreducibility of C, there are two elements %q G M and jo £ M C such that 
Vi -+j > 0; remembering that, by definition, if % ^ j then Vi-*j < £i-^j, we have, 



> 



1 ^ ] ^io-^j 
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JO 



*o->Jo 
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> ^ f 1 - 5°. ) >0 . (7.42) 



And this proves the regularity of P; from the Markov theorem follows the uniqueness 
of the invariant probability distribution and this proves the theorem. □ 



In most practical cases, Eq. (7.40) for i ^ j is written in the form 



7^ = £^- mini 1,^1 . (7.43) 



7T,; 



The meaning of this relation is that the entire Markov chain can be built with pre- 
determined moves, that might be accepted with probability min ^1, These 

moves, starting from a probability X n propose a transition to a probability X„ +1 : if 
this transition is accepted, X n+1 has been determined, in case of rejection X n+1 = X n . 
The condition of irreducibility of £ means that the moves must be chosen so that their 
combination is able to explore the whole state space E, this property is called ergod- 
icity. The symmetry of £ is a detailed balance condition on the Metropolis moves, 
this however can be dropped in favor of the weaker condition C-^j > whenever 
Cj^i > if a new definition of is taken: 



V, C, ..minfl.^ A . (7.11) 

A good choice of Metropolis moves will enhance the convergence of the Markov chain 
towards the equilibrium probability distribution. 
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7.1.3 Sampling and expectation values 

In the previous section, we have seen a procedure to sample an arbitrary distribution 



probability such that in Eq. (7.3). In order to evaluate that iV-dimensional integral, 
however, there are still two problems to take care of. First, the Monte Carlo evalu- 
ation of an integral is a statistical method and as such the results have an intrinsic 
statistical error due to the finite number of sampled values, and second, the Metropo- 
lis algorithm produces an highly correlated sampling of the probability distribution 
p(x). To overcome these problems the block average technique can be used: if the 
summation in Eq. (7.1) is truncated after N s terms, a block In s = jj~J2f=i 9(%i) 



is defined. If N s is long enough, each evaluated In s can be considered statistical 
independent from the others; hence, by the central limit theorem, is a Gaussian 
distributed random variable; if the set of values 1^,1^, I^ blocka are generated with 



Eq. (7.1), then the average value I avg and the standard deviation a can be computed 
with the usual formulas: 



Nbloc. 



a 2 = - 

Nblocks 

i=l 



-t * ' OLOCKS 

i avg = — Yl p n ( 7 - 45 ) 

blocks ■ , 
i=l 

I ^blocks 

— 7 E ( J n-w> 2 > ( 7 - 46 ) 



where the error of the estimation I avg is a /y/Nuocks- This method for the evaluation 
of TV-dimensional integrals is asymptotic, in fact the correct sampling of a distribu- 
tion will be given only after that a long enough Markov chain has been built. This 
means that the probability distribution p(x) is sampled only after a certain number 
of equilibration steps. This equilibration number can be evaluated by plotting on a 
graph the value of the integral averaged within a Monte Carlo block versus the index 
of the corresponding Monte Carlo block: equilibration is over as soon as transients 
disappear from the plot, provided that the chosen set of moves can efficiently explore 
the whole space of events. This might seem to be an easy condition to fulfill but in 
some cases a long equilibration is required, especially when the distribution density 
to sample has many local maxima separated by regions of low probability density. In 



Fig. 7/2 we show the value of an integrad evaluated at each MC step: the equilibration 
transient is clearly visible in the first twenty MC steps; the correlation of the Markov 
chain manifests here as a pattern in the values of the integrand. 

7.1.4 Metropolis sampling 

We now apply the Monte Carlo sampling to the proble m of evaluating a quantum ex- 



pectation value of a local operator O introduced in Sec. 2.1 this quantum expectation 



value, from Eq. (2.11), can be written in compact form as 

'd\ = [ dTO(T)p(T) (7.47) 
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Figure 7.2: The "istantaneous " value of an integrand for each Monte Carlo Step. 
This integral represents the energy of one atom of 4 He in a ID model potential defined 
by V(x) = a\x A — (T2X 2 , where <j\ = 8 KA~ 4 and (j 2 = 8 KAr 2 . The methodology 
used to evaluate this integral is the Path Integral Ground State that will be described 
in Sec. I2.ll 




where T = {R\, Rm}- In this section we consider the PIGS case; the adaptation to 
quantum thermal averages is straightforward and will be considered contextually. The 
Metropolis algorithm is used to sample the multi-dimensional probability distribution 



p(T) in Eq. (2.18), that, explicating the normalization constant A/", takes the form 



frr (gi) YlfJi 1 G (Rj, R j+1 , St) (Rm) 
j dT ^ T (R ) n'lT 1 G (Rj, R j+ i, r) * T (R M ) 



P(T) 



(7.48) 



The sampling of p(T) is made with a sequence of Metropolis "moves". A move is 



a two-step process sketched in Fig. |7.3| this process, from a set of configurations T 
proposes a new set T new and then evaluates whether to accept or reject the new set 



of configurations. The probability to accept the move is defined by Eq. (7.43), where 
the term Ci->j represents the probability to try a move that from a configuration % 
proposes a new configuration j. From the a-priori knowledge of the system under 
study it is possible to use guided moves that are more likely able to sample physical 
configurations rather than highly improbable ones (in this case the probability to 
accept the move becomes Eq. (7.44); such a guided approach would enhance the 



convergence of the sampling, especially if the probability distribution has many local 
minima. The moves that will be described shortly are unguided, so that C^j = Cj->i 
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Figure 7.3: General scheme for a Metropolis move: from a configuration (a), a move 
is proposed (b). At this point the move can be accepted (c) or rejected (d). If the 
move is accepted, the new MC step will have a new configuration T new ; otherwise the 
same configuration r o ^ is sampled again. Grey beads and lines represent the removed 
segment of the polymer. 



and their probability to be accepted simplifies to the following relation 



a \Lnew) = mm 1, — . 7.49 

^TmXlfj'GiR^R^dT^TiRM) ) 

In the context of quantum-classical isomorphism, a move can involve one or more 
different polymers. A move involving a single polymer is represented as a reconfigu- 
ration of some or all the beads of the polymer itself; the probability to accept such a 
move depends on the correlations of the beads at the same imaginary-time discretiza- 
tion (inter-polymer correlations) and on the correlations between adjacent beads that 
belong to the same polymer (intra-polymer correlation). This latter contribution is 
also referred as kinetic spring because comes from the kinetic term of the Hamiltonian 



(Eq. |2.13 ) and disfavors configurations in which two adjacent beads are placed far 
away from each other. A move involving many polymers is a generalization of a single 
polymer move; these moves can also involve permutations between polymers that can 
be employed to take into account the quantum statistics of the system. 

A Monte Carlo simulation consists of a set of Monte Carlo Steps (MCS); in general, 
after each step the estimators can be evaluated. A MCS consists of a set of Metropolis 
moves that are tuned so that the effect of all the accepted moves modifies the positions 
of roughly half the beads that compose the system of polymers. The Metropolis moves 
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that are proposed here are the translation moves, the Brownian bridges and the 
permutation sampling. Later in this section will be introduced the Worm algorithm 
in the Canonical ensemble with its Metropolis moves. These moves can be also used 
for PIMC with a slight adaptation for the translation of a polymer. 

In order to simplify the notation, we define the free particle propagator that 



appears in Eq. (2.13) as 



Go 



1 



(7.50) 



Translation of a single bead 



In a PIGS simulation, this move, represented in Fig. 7.4 is generally applied to the 
first or the last bead of a polymer, namely r\ or rf 1 . This is the simplest move 
for such beads: there are other, more performant, possibilities that allow to move a 
certain number of beads including f\ or rf f ; one of those moves is a generalization 
of the Brownian bridge that takes into account the correlations from the trial wave 
function. The Brownian bridge will be introduced soon; however, in this work we 
did not implement the mentioned extension. Instead, we used the following move. 
Let's focus on the first bead of the z-th polymer, r\\ the new configuration will have 



a translation of a vector d applied to f] 



namely r] 



fj + d. 




Figure 7.4: Scheme for the translation of an extremal bead of a polymer. Grey beads 
and lines represent the old position and kinetic correlation of the bead. 



The probability to accept this move is 



sing 



a({R\ 



min (1, Pi 



j 1 tr) 



-££ Mi <|r4 ew -rt|) 



2 



(7.51) 
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For the last bead, ff , 



tr) 



Psing — St ^ /i „ • ( 7 -52) 



e ^ Efe ^KK W ^ \)ty T (R°M) 
Translation of a polymer 



In a translation (Fig. 7.5 a polymer z is rigidly moved by a vector. Kinetic springs 
remain unchanged and the probability to accept the move depends thus only on the 
interpolymer correlations. 




Figure 7.5: Scheme for the translation of an entire polymer. Grey beads and lines 
represent the old configuration of the polymer. 



a({R} r 



min (1, Pi 



2 T T.k^i "Genera r k\ ) 



tr 



m T {R°i d )e 



tr) 



x 



X 



Hf = - 1 e - Sr £m I ) e -£ «( | n£» I ) ^ T (i^) 



(7.53) 



yold\ 



In PIMC, the translation move has a probability to be accepted that is slightly 



different from Eq. (7.53): 



a({R} new )=mm(l,P tr ) 



p 



1) 



tr 



-St y~\ / - v(\r 3 



(7.54) 



125 



Brownian bridge 



The Brownian bridge move is a very efficient way to reconstruct a segment composed 



of S adjacent beads. A schematic description of this move is shown in Fig. 7.6 
The segment of polymer that the Brownian bridge re-creates represents the sampling 
of the free particle propagation in imaginary time between two time sectors. This 
propagation, as will be shown below, can be sampled exactly via the Box-Muller 
method[0]. This is a useful feature: the kinetic correlations of the reconstructed 
segment are sampled exactly, therefore the probability of acceptation of the move 
will depend only on the correlations between different polymers. This is readily seen 
if one considers the correlations of a segment of the i— th polymer (r^ 7 , ...,r? ) that 
is part of a configuration of polymers V: 

TT(r) = f[ e -^M m ^T +1 \\-^M\rT-rr\) (7.55) 
m=j 

the kinetic and potential parts are factorized. Recalling the Metropolis algorithm 



(Sec. 7.1.2), the general form of the transition matrix that represents the Markov 
chain is factorized in a "move" and an acceptance of the move: 

V^^C^mmfl, 7 ^^) . (7.56) 




Figure 7.6: (Upper panel) The Brownian bridge: a segment of a polymer is recon- 
structed. (Lower panel) Iterative procedure used to sample a free particle propagation 
between two fixed extremities. Grey beads and lines represent the removed segment 
of the polymer. 
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In the case of the Brownian bridge, the move Ci-+j = Cj-n is the exact sampling 
of the free particle propagation in imaginary-time, Hj = 7t(T new ) is the correlation 
value of the new segment and 7Tj = 7r(T ^) is the correlation value of the segment 
before the move: the rest of the system is not changed by this move and the respective 



correlations, being unchanged, cancel out. From Eq. (7.55) it is clear that the free 
particle propagator that defines C^j is a part of ttj; therefore, the probability to 
accept a new segment becomes a = min (1, Pbb), with: 

P bb = ^ 7j K (7-57) 

1 Lm=j 

where the reconstruction starts from the j + 1 bead of the i— th polymer and the last 
reconstructed bead is at position j + s — 1. 



The following operations, also illustrated in Fig. |7.6[ are performed during the 
move: 

• Remove the beads between j and j + s. 

• Create a new timeslice at position j + 1: from the coordinates of the timeslices 
j and j + s, determine the coordinates of the new timeslice j + 1. These 
coordinates are determined in the following way: we first note that the 
bead at position j + 1 is the free propagation from the bead at position j, 
Pi(r?, r*) = Go(rf , f*, 8r), and the free propagation from the bead at position 
j + s, P2{?i +i \ r*) = Go(r? +s \ r*, s5r). The probability density from which the 
position r? +1 is sampled is thus the joint probability p±p2', with straightforward 
algebraic operations, this joint probability can be reconduced to the Gaussian 



form of Eq. (7.50) times a trivial normalization constant Mt that won't affect 



the sampling; namely, the new bead is sampled from 

p{?*) = Pi(r7,r> 2 (r7 +S ,f*) = Af t G (r\f? + ^ ~ r * , j^Sr) (7.58) 



• From the newly created bead, and the bead j+s, determine the coordinates 
of the bead j + 2. This is done by iteration, considering a segment of polymer 
that starts at j + 1 and has length s — 1. The procedure is iterated until the 
bead at position j + s — 1 has been determined. 

This move samples the free particle propagation between two given extremities; 
it may happens that the two extremities are in different polymers that are connected 
at a timeslice j p by a permutation cycle P; in this case the labels j, i and k must 
be permuted, so that % — > Pji and k — > Pjk, where Pj = I for j < j p and Pj = P 
otherwise. 

The probability of acceptation can be varied modifying the length of the Brownian 
bridges. As an empirical rule, a good choice of this probability can be between 0.3 
and 0.5. This is a reasonable trade-off between long moves and small moves: long 
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moves, on one hand, would yield an high rejection rate resulting in poor performance 
of the simulation; small moves, on the other hand, have an high acceptance ratio but 
might compromise the ergodicity of the simulation; this happens because unprobable 
configurations would rarely be sampled. 



Permutation sampling 

The permutation sampling introduces the Bose symmetry in the system. In the 
polymer description, it is a move that involves a permutation between a variable 



number of polymers greater than one. As mentioned in Sec. 2.1.1 this move is 
necessary in PIGS when the trial wave function does not possess the Bose symmetry. 
This move is also necessary at finite temperature: in PIMC, in fact, the thermal 
average is not expressed through the quantum imaginary-time evolution of a trial 
wave function; as consequence, the Bose symmetry has to be explicitly introduced 
through permutation sampling. 

Here we will show the algorithm described in Ref. [TU]. The polymers involved 
in the permutation are selected with a kinetic test to be described soon. Once the 
polymers i\ and i<i have been selected, their beads between two time sectors jo and 
jo + s are removed. At this point, a Brownian bridge is made from the bead at position 
jo of the polymer z 1; to the bead at position jo + s of the polymer i 2 . The permutation 



move follows the scheme in Figure 7.7 the iteration of this 'swap' procedure proceeds 
between polymers %2 and i^, and so on until a polymer i n closes on the polymer i . 




Figure 7.7: Permutation of two polymers i and j: the resulting configuration has 
still two polymers of the same length; it is thus topologically similar to the previous 
one. Grey beads and lines represent the removed segment of the polymer. 



There are two main steps in permutation sampling: the kinetic test and the re- 
construction step. 

Kinetic Test step. Given a starting timeslices jo an d a length of reconstructions 
s, this operation selects the polymers that are best suited for permutations and gives 
in output an ordered sequence of swaps between the polymers. The kinetic test starts 
from a random polymer i\ that is chosen by generating an integer number between 1 
and TV from a uniform distribution probability. At this point, the following operations 
determine the next polymer which joins the permutation cycle. 
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For the particle ii, build a table as follows 

Kl» = °o (ft ^ s&r) (1 - 6 h J (7.59) 
where r?° are the coordinates of the polymer i\ at timeslice jo- 



From Eq. (7.59), the probability to accept the particle i\ in the permutation 
cycle is 

C (1) = —. T . (7.60) 

Generate a random number p, uniformly distributed between and 1. If p > 
the move is rejected; otherwise the process continues to the next step. 

From K\ u a new particle is randomly chosen. The probability to chose a the 
particle v is 

n„ = =%- . (7.6i) 

The new particle is selected with a 'faked roulette': the interval [0, 1) is parti- 
tioned with bins of width U v ; a bin U u corresponds to the particle v\ generate 
an uniformly distributed random number in the interval [0,1), the bin which 
contains this number corresponds to the particle that 'wins' the faked roulette; 
this particle is 

Make an acceptance test on i 2} similarly to that made on %\ 



KL = G (fl ,r^ +s ,s8r) (1-8^) 
Go ^,r^ +s ,s5r)+^Kl 



C (2) = —. T . (7.62) 



Again, generate a random number p, uniformly distributed between and 1. If 
p > C^ 2 \ the move is rejected, else the particle i 2 is added to the permutation 
cycle. 

These operations are repeated until either an acceptation test fails or a particle i a = i\ 
is added to the permutation cycle. 

The Dirac's deltas that appears in K™ nU have two purposes: 

• Exclude the possibility for a particle to swap with itself 

• Exclude the possibility for a particle to swap with any other particle already 
added to the permutation cycle, except for the first particle of the permutation 
cycle. 
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The general definition for Kf , , is thus 

K« = Go (t£, e +s , s5t) (1 - 5 l2 , w ) (1 - 8 i3 , w ) ... (1 - 5^) . (7.63) 

The output of the kinetic test step is a sequence of particles (zi, 12, i a , h)', from 
this output, the reconstruction step begins. 



Reconstruction step. Starting from the previously obtained permutation cycle 
(ii,i2, ■■-,i a ,ia+i = h), ol Brownian bridges are built from the j bead of the i\ poly- 
mer to the jo + s bead of the i%, an so on until the last Brownian bridge from the jo 
bead of the i a polymer to the jo + s bead of the i a +i closes the loop. 

The obtained new configuration has a probability to be accepted that is the prod- 



uct of Eq. (7.57) for each reconstructed segment. If the move is accepted, this new 
configuration is kept, if this acceptance test fails the configuration prior to the permu- 
tation move has to be restored. The probability to accept exchanges is usually very 
low and in order to obtain an efficient permutation sampling one has to try thousands 
of permutation moves in a single MCS; moreover, in most cases, the probability to 
accept a permutation drops exponentially with the number of polymers involved in 
the permutation and thus the efficiency of this algorithm for permutation sampling 
get worse with increasing particle number N. Given a permutation cycle, the prob- 
ability to accept the reconstruction step is roughly the product of the probability to 
accept each single Brownian bridge of the same length; this suggests that in most 
cases, a good choice of the length of reconstructions s can be roughly the same as 
that of a single Brownian bridge; this however may not be true if the polymers are 
distant each other; in this case the probability to accept a permutation is maximized 
if it involves a sufficiently large imaginary-time; this holds even though the Brownian 
bridges in the reconstruction step would have low acceptances. 

7.1.5 Estimators 

Let's consider again the expectation value of a local operator O 

o) = J dTO(T k )p(T) (7.64) 

where T = {Ri, Rm } and 1 < k < M represents the position in the path integral 
at which the operator is applied. This equation holds for both PIGS and PIMC 
depending on the choice of the multi-dimensional probability distribution p(T), to be 
more specific, in the PIMC case, 

P(T) - P™ C (D = J|g (7.65) 
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where we set Rm+i = Ri- In the PIGS case, 



^t(Ri) nJlT 1 G (Rj, R j+ i, St) V t (R m 



p(T)=p PIGS (T) 



M ) 



J dT * T (i20 n'lT G (R^Rj+i, St) Vt(Rm) 



(7.66) 



In the previous Sec we described a method to sample p PIGS (T) and p PIMC (T), here 
we focus on the application of the operator O to the density matrix. We have al- 



ready pointed out in Sec. |2.2| that in PIMC, due to the cyclic property of the trace 
operation, one can shift the position of O along the path integral without chang- 
ing the expectation value. This is useful because one can use all the configurations 
{Ri, Rm} to compute the expectation values and then average the results. Also 
in PIGS an operator can be evaluated at any imaginary-time 77; to be more explicit 
on the meaning of "application of an operator at an imaginary time 77": 



(7.67) 



q? ( T = ISt) |0|* (t = (M - 1)St) 
J dRi...dR M ®t (Ri) ...G(Ri-i, R h St)6 (R t ) G{R h R l+h St)...V t (R m ) 



J dR 1 ...dR M #t (Ri) G{R 1 ,R 2 , St)...G{R M -i, Rm, St)^ t (R 



M 



with 2 < I < M — 1. Here |^(r)) represents the evolution of the trial wave function 
|^t) at a n imaginary-time r, namely \^(t)) = \e~ TH ^T)- Differently from PIMC, 



due to Eq. (2.7), only for r < 77 < r — r it is verified that, to a good approximation, 
|tf(r = n = ISt)) ~ |0) and |*(r = r M _/ = (M — 1)St)) ~ |0); in this case, Eq. (T67\) 
becomes an expectation value on the ground state of the system. Outside the interval 
[r ; r — r ] the expectation values are mixed, more specifically, for an imaginary-time 



index h so that T h < r , (\I>(r = 



O 







For 7/j = and Th 



> and ( 



O 



t we obtain 
with the trial 



respectively the mixed expectation values y& t 
wave function ty T . These expectation values are obtained by applying the operator 
directly on the trial wave function, 



/ dR x ...dR M 6 (i?i) #r (Ri) G(R u R 2 , St) 



.G(R M -i,Rm,St)V t (Rm) 



(7.68) 



/ dRi...dR M (Ri) G(Ri,R 2 , St)...G(R M -i, Rm, St)^ t (R 



M 



an analogous relation holds for ^ (Rm)- The mixed expectation value is very use- 
ful in the evaluation of the total energy. The Hamiltonian H , commutes with the 



propagator e SrH and thus the total energy can be evaluated at any time-step; how- 



ever, in Eq. (7.67), we are approximating e with a propagator G(R, R' , St); as 



consequence, the commutation rule that allows the evaluation of the total energy at 
any time-step, holds only if the small-time approximation of the propagator is accu- 
rate enough; this provides a useful check of convergence of PIGS for what concerns 
the choice of St. Moreover, if an accurate trial wave function is available for the 
system to study, the evaluation of the total energy on the variational wave function 
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is more accurate than that at other imaginary times; this happens because the trial 
wave function introduces correlations that guide the Metropolis sampling in a more 
efficient way to explore physical configurations. 

Now that the effect of the imaginary time evolution on the expectation values 
has been described, we show the application of the operator O on a density matrix 
G (Rj,Rj + i, Sr) in the most common choices of O. This, of course, is generally depen- 
dent on the particular choice of small imaginary-time approximation for G, exception 
made of the case of an operator that is diagonal in the coordinate representation, for 
example, the static structure factor, the density-density correlation function in imag- 
inary time or the radial distribution function. Consider the general case for a density 
matrix that can be expressed as follows: 

G (Rm, Rm + i,Sr) = G (R m , R m+1} Sr) e-^-ViA) ( 7 . 69 ) 

where Go is the density matrix for free particles. Some of these density matrices are 
illustrated in Appendix. [Bj Let's also consider, for simplicity, the PIMC case, so that 
the partition function of the system becomes 

/M _ 2 

n dR m e- ^ +l) e SrU(R m , Rm+1 ,Sr) . (7 70) 

m=l 

The following discussion applies also in the PIGS case with at worse slight modifica- 
tions that will be described contextually. 



Energy 

The Energy per particle, E/N, is the expectation value of O = H/N. Apart from 
the physical importance of this quantity, the energy is one of the main expectation 
values that are used to tune the parameters of the simulations. 



Hamiltonian estimators This estimator is obtained by applying the operator 
H = T + V to the density matrix (7.69). 

The potential term V = J2i<j v ( r *i) i s diagonal on coordinates representation, so 
that 



VG (R m , R m+1 , Sr) = V (R m ) G (R m , R m+1 , Sr) 



(7.71) 



In Fig. 7^ we show the potential energy obtained from a PIGS simulation of two- 
dimensional 4 He. 

The kinetic term is T = — ^- ^1 an d 



V 4 2 G (R m , -R m +i) St) — G (R m , -R m +i) Sr) 



(7.72) 
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Figure 7.8: PIGS computation of the potential energy per particle versus imaginary- 
time for a 2D system of N = 16 atoms of 4 He at a density p = 0.045 A~ 2 , interacting 
with the Aziz potential described in Ref. [7]. The trial wave function is >— 1 and 
the total projection time r = 0.5 K -1 . Red squares are obtained with an 8-th order 
multi-product expansion (see Appendix [b|) at a timestep 65t = 1/80 K -1 . Black 
circles are obtained with the primitive approximation at a timestep 5t = 1/480 K -1 ; 
however, for comparison purposes, only points at r m = 6m jbr are shown. 



where U = U (R m , R m +i, St) for simplicity, A = |- and d is the dimensionality of the 
system 



Thermodynamic estimators The total energy per particle can be obtained also 
from the thermodynamic definition, 

E(N,V,P) 1 dZ(N,V,P) 

N ~ NZ d(3 {L<6) 



where Z is the partition function defined in Eq. (7.70). The thermodynamic estimator 
for the energy per particle is thus 

A / </ ... x - ,-. - ■ •." >.'•/))• »•.»+!••" i \ (774) 



N \ 2St AXSr 



1 V^^p 7~> x2 I dU (R m ,R m+1 ,5T) \ 

■2 MN 2^^ n ™ Km ^> + MN d5r / 

m=l I 
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In the same way, the kinetic energy per particle, K/N, can be obtained from the 
thermodynamic relation 



K _ m dZ (N, V, (5) 
N ~ JZ dm" 



(7.75) 



and the estimator becomes 



K 

N 



d 



M 



1 £ {Rm - Rm+1 f + s Zj U(Rm l Hm+lM) > (7-76) 



25t AX5r 2 MN^ y "" 1 ~~ I " rx ' '5tMN dm 

m=l 



where (...) is the average on the configurations {R m }^ =1 that are sampled by the 



metropolis algorithm. An example of application of this estimator is shown in Fig. 7.9 

l 




Figure 7.9: PIGS computation of the kinetic energy per particle versus imaginary- 
time for a 2D system of iV = 16 atoms of 4 He at a density p = 0.045 A -2 , interacting 
with the Aziz potential described in Ref. [7]. The trial wave function is >= 1 and 
the total projection time r = 0.5 K _1 . Red squares are obtained with an 8-th order 
multi-product expansion (see Appendix [b|) at a timestep 65t = 1/80 K -1 . Black 
circles are obtained with the primitive approximation at a timestep 5t = 1/480 K _1 ; 
however, for comparison purposes, only points at r m = 6m /St are shown. 



The Hamiltonian and the thermodynamic estimators provide two different ways to 
obtain the energy, however they suffer from statistical fluctuations that increase with 
smaller values of St, this is particularly true for the Hamiltonian estimator due to 
the presence of the laplacian operator, but happens in smaller degree also in the 
thermodynamic estimator because the first two terms in Eq. (7.74) and (7.76) are 
quantities that increase when St decreases and cancel each others. This requires 
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longer simulations when St is small and usually poses a computational limit for the 
evaluation of the total energy. There are at least two possibilities to overcome this 
problem: one can either use higher order estimators which achieve convergence at 
higher timestep or introduce a more advanced estimator. A choice for the latter 
possibility is the virial estimator [TT]; the derivation of this estimator is shown in 
appendix |A| together with an explicit derivation of the thermodynamic estimators for 
the Pair Suzuki approximation. In this appendix we show also that although the 
thermodynamic estimators are obtained from a thermodynamic relation, because of 
the similar formalism of PIGS and PIMC, it is possible to use these estimators also 
in PIGS. 



Radial distribution function 



The pair correlation function g (r*i,r*2) is the probability to have a particle at f\ and 
a particle at r 2 . Within the path integral formalism, 



9 \r\,T2) 



z 



dr 3 ...dr N G(R,R,fi) 



(7.77) 



In a uniform system the pair distribution function depends only on the distance 
r = | r\ — r*2 1 , with a change of integration variables in Eq. (7.77) and using the 



definition of thermal average (7.64), the estimator becomes 



9{r) 




(7.78) 



where we have taken into account the symmetry under particle exchange and the esti- 
mator has been averaged over the timeslices m in order to employ larger statistics. In 
PIMC the sum over m covers all the timeslices; in PIGS, this sum must be intended 



only over the central timeslices, where Eq. (7.67) gives an accurate description of the 



ground state. To evaluate this estimator in a computer simulation one defines a par- 
tition P n of the interval [0; Li/2] where Li is the smallest size of the simulation box, 
and every element of the partition has a length Ar; with P n = [nAr; (n+ l)Ar] then, 
construct an histogram of the frequencies of the relative distance between two 
particles of the system at the same imaginary-time index. This histogram has to be 
normalized with the number of particles that a free particles system of the same den- 
sity would have at a bin n, namely for d = 3, V n = y^n [({n + 1) Ar) 3 — (nAr) 3 ] . An 
example of QMC evaluation of the radial distribution function for a two-dimensional 



system of He is provided in Fig. 7.10 



Static structure factor 

The static structure factor is useful to study the spatial order of a system in the 
reciprocal lattice; it is in fact connected to the pair distribution function by a Fourier 
transform. This estimator is defined as a quantum average of the density operator 
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Figure 7.10: PIGS computation of the radial distribution function of a 2D system of 
N = 16 atoms of 4 He at a density p = 0.045 A -2 , interacting with the Aziz potential 
described in Ref. [7j. The trial wave function is >= 1, the total projection time 
t = 0.5 Kr 1 and the averages were taken in the imaginary-time interval 0.2 K _1 
- 0.3 K _1 . Red squares are obtained with an 8-th order multi-product expansion 
(see Appendix [b]) at a timestep Q5t = 1/80 K _1 . Black circles are obtained with 
the primitive approximation at a timestep 5r = 1/480 K™ 1 . The simulation box is a 
square of late L; the radial distribution function has been computed also in the range 
{L/2-LV2). 



s k 



N \ p * p -*/ ~ NZ 



N 



N 



dR p (R, R,/3)[J2 ^ E ^ ■ ( 7 ' 79 ) 



,i=i 



. i=l 



Using the Euler identities, the static structure factor can be expressed in a form that 
is computable by different Monte Carlo methods 



S ( k 



1 



N M 



NM 



2^ cos ( k ■ f™ \ cos ( k ■ Tj m j + sin ( k ■ fl m ) sin ( k ■ r 



izfzj m=l 



where we have averaged over the M equivalent timeslices as before. As in the 
previous case, in PIMC the sum over m covers all the timeslices; in PIGS, this sum 
must be intended only over the central timeslices, where Eq. (7.67) gives an accurate 
description of the ground state. 
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Figure 7.11: PIGS computation of the static structure factor along the axis direction 
of the simulation box. The system is two-dimensional and consists of N = 16 atoms 
of 4 He at a density p = 0.045 A~ 2 , interacting with the Aziz potential described in 
Ref. [?]. The trial wave function is >= 1, the total projection time r = 0.5 
K _1 and the averages were taken in the imaginary- time interval 0.2 K" 1 - 0.3 K -1 . 
Red squares are obtained with an 8-th order multi-product expansion (see Appendix 
[b]) at a timestep QSt = 1/80 K _1 . Black circles are obtained with the primitive 
approximation at a timestep St = 1/480 K _1 . 



It must be remarked that in a QMC simulation the simulation box has finite 
dimensions that in the case of d dimensions are (L , L^). This implies that the wave 

vectors that are accessible by the simulation are of the form k = ^f^o ; jr d Ud 



with n , ...,n d G I. An example of static structure factor is shown in Fig. 7.11 



Imaginary— time correlation functions 

The static structure factor is also called 'density-density correlation function'. In 
general, a correlation function between two operators A and B is a quantum average 

c AB = (Ab} . (7.80) 

In particular, the operators A and B can be evaluated at different imaginary-times 
(different time-sectors) and this is an imaginary-time correlation function. A signi- 
ficative example of imaginary-time correlation function is the density-density one: 

F(£,T)=i(ft(0)p_ 5 (T)) (7.81) 
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that for r = reduces to the static structure factor. This function is related to the 
dynamical structure factor of the system by a Laplace transform and thus contains 
informations about the excitations of the system. However, these informations are 
accessible from Eq. (7.81 ) only by solving a numerical inverse Laplace transform; this 
operation is carried out on a function F(k, r) which is known only for some imaginary- 
time r m and with a statistical uncertainty. In these conditions, the inversion of the 
Laplace transform is an ill-posed problem. Methodologies to face those problems 
have been implemented by many research groups [ 



One body density matrix 

An important quantity to sample is the one-body density matrix (OBDM) because 
it is strictly connected to the Bose Einstein condensation (BEC). 

BEC can be defined as a macroscopic occupancy of a given quantum state. The 
simple BEC occurring in 4 He may be described by a momentum distribution of the 
form 



n(p) = N 5 {p) + n(p) 
The OBDM of the ground state |0) of the system is 



Pi[r,r 











(7.82) 



(7.83) 



The Fourier transform of this equation gives the momentum distribution of the 
system at its ground state: 







alap 







a ,7 



dre* p ' r ^ (r) 



(7.84) 
(7.85) 



Placing Eq. (7.85) in Eq. (7.84) yields 



n. r . 



drdr e 



■ip-ir-r 



(2nhy 
i 



, i drdr e hV v r > p\ [r,r 



{2-nhy- 



S -> s 

, . dtdse~^ ps pi ( t + -,t — - 



(7.86) 



where s — f—f and t = (r + r ) /2. If the system is uniform and isotropic, the 
OBDM depends only on s — and in the thermodynamic limit 



n,- 



V 
(2irhf 



dse '^' g pi (s) 



(7.87) 
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If the momentum distribution in Eq. (7.82) is anti-transformed to coordinates, an 



OBDM with unvanishing tail at high s is obtained. An OBDM that displays such 
asymptotic behavior, that is lim^oo pi (s) = n > 0, has an off-diagonal long range 
order; if that function is normalized so that p\ (0) = 1, n corresponds to the fraction 
of BEC. 

The OBDM in the path integral notation is 

pf /GS (r,f') = J dR l dR M df 2 M/2 ...df^ 2 ^ T {R 1 )G{R l ,R M/2 ,r/2) x 

xG(R m/2 ,Rm,t/2)^ t (R m ) (7.88) 

where R m = (f, r 2 m , r™), R m = (f ', r 2 m , f ™) and dR m = YliLi dr™. In the 
PIMC case, an analogous procedure gives 

p piMc ^ /j = Vfdr 2 ..df N G(R, R',P) ^ (7 gg) 

In the classical isomorphism, the z-th polymer has been split: if the original polymer 
had M beads, the bead at a position j is removed and two new beads are inserted. 
In PIGS j should correspond to an imaginary-time projection large enough to have 
an accurate description of the ground state; in PIMC j can be anywhere in the path 
integral. These beads, f and r , are not linked each other but are respectively the 
last and the first timeslice of the new open polymer that would appear. In the PIMC 
case, this will be a single open polymer with extremities r and r ; in the PIGS case 
two half polymers will appear, with extremities fl,f and r ,r*j M . 

If the system is homogeneous, p\ depends only on the distance r = \r — r\; the 
OBDM is then sampled making an histogram of the relative distance that the two 
newly created beads have at every MCS. This histogram should be normalized so 
that pi(r = 0) = 1; this is done after the QMC evaluation of pi with a Gaussian 
fit of the small r part of pi{r). In some cases, this might not be the optimal choice 
for the normalization: the Worm algorithm, described in Sec. |7.1.6 can provide an 



a priori normalization of the OBDM with a QMC evaluation of M{Z). As example 



we show in Fig. 7.12 the one body density matrix for a two-dimensional system of 
4 He. We point out that in the PIMC case, the sampling of permutations is crucial 
to obtain off diagonal long range order: permuting ring polymers with the open 
polymer will result in a bigger open polymer that allows its extremities to get far 
away each other, eventually contributing to a non vanishing tail in the OBDM. In 
PIGS, permutations will yield simply other open polymers and in most situations 
are not essential; however, there are cases in which permutations are essential in 
order to obtain an ergodic sampling of the configurations; we have indeed shown an 



example in Sec. 2.1.1 when we studied, with PIGS, the condensate fraction of 4 He 
with a particular choice of ^t'- a Gaussian wave function centered on the equilibrium 
positions of solid HCP 4 He; this is a very particular case in which the trial wave 
function introduces correlations that constrain the terminal beads of the polymer 
to arbitrary positions. In this context, single polymer moves will have the same 
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Figure 7.12: PIGS computation of the one body density matrix of a 2D system of 
N = 16 atoms of 4 He at a density p = 0.045 A -2 , interacting with the Aziz potential 
described in Ref. [7]. The trial wave function is >— 1 and the total projection time 
t = 0.5 K _1 . Red squares are obtained with an 8-th order multi-product expansion 
(see Appendix [b]) at a timestep 65r = 1/80 K -1 . Black circles are obtained with 
the primitive approximation at a timestep 5t = 1/480 K _1 . The simulation box is 
a square of late L; the radial distribution function has been computed also in the 
range (L/2;La/2). The normalization constant has been computed with the Worm 
algorithm (see Sec. 7.1.6). 



pathology of the PIMC case: the extremities of the two half polymers are pinned, as 
consequence a move that would increase the distance between r and r will be soon 
rejected by the stretching of the kinematic correlations; a permutation with another 
polymer, on the other hand, will allow the sampling of long-range order exactly as in 
the PIMC case. 

In some situations, the OBDM decays exponentially; this is typical in most solid 
systems. In such contexts, a common way used to sample the OBDM at large dis- 



tances r is with the introduction of a repulsive factor f(r) in Eq. (7.88) (or Eq. 7.89 
for PIMC) 



f(r) 



1 



1 + Ae~ Br2 + Ce~ Dr 



(7.90) 



The parameters B and D are tuned to fit the exponential decay of the OBDM whereas 
A and C determine the strength of the repulsive factor. The repulsive factor is aimed 



to modify the probability densities (7.88) and (7.89) to a roughly uniform probability 



density p(r) that is easier to sample. Once the histogram of p(r) has been obtained, 
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the OBDM p(r) is recovered with a reweighting of p(r): each histogram bin (p(rj),rj) 
is divided by the quantity /(r^). 



Superfluidity 

It has been derived in Ref. pQ that the superfluid density of a system can be expressed 
through the winding number W, 



(7-91) 
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" 2A/3iV 


and the winding number 


is defined by 






1 N 

i=i 


'o 


dfi(ty 

dt 



(7.92) 



In the path integral notation, the discretized expression for Eq. (7.92) becomes 

M N 
m=l i=l 



w 



m+l 



(7.93) 



where L is the late of the simulation box, r^ I+1 = . and P\ the permutation 
operator introduced for the sampling of the Bose symmetry. The winding number 
represents the number of polymers that wind around the periodic boundaries of the 
simulation box. Averaging over the d spatial dimensions, the winding number esti- 
mator in periodic boundaries conditions is 



Ps 

P 



W' 



2d\(3N 



(7.94) 



Superfluid density at T = K At zero temperature the superfluid fraction can 



be obtained with the center of mass diffusion in imaginary-time^^. Eq. (7.94) in fact 
can be viewed as the ratio between the diffusion constant D c of the center of mass of 
the system and the diffusion constant of the non-interacting gas, Do = h 2 /2m. The 
diffusion of the center of mass is obtained from the long r limit of this relation 



lim — 

T x ^OO 4 



Rcm (r,) - R CM (0) 



(7.95) 



where the center of mass at a discrete imaginary time r 
Eli rT/N. 



mdr is Rcuij) 
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The superfluid density becomes 



Ps 

p 



Dc 

Dn 



lim — 

T— >0O 4A 



R CM (r) - Rcm (0) 



(7.96) 



In a PIGS calculation one has to make some considerations: first, this estima- 
tor can be applied only in an interval of imaginary-time [r , r — r ] so that, in or- 



der to achieve long-r convergence in Eq. (7.96) one has to employ sufficiently long 



imaginary-time projections; secondarily, the PIGS method does not explicitly fix the 



center of mass of the system and thus Eq. (7.96) cannot be used when the center 



of mass of the system is allowed to drift; this happens, for instance, in bulk homo- 
geneous systems, where the system can translate freely in the simulation box: due 
to the property of PIGS to give an unbiased sampling, any correlations from the 
trial wave function that would eventually fix the center of mass is removed by the 
quantum imaginary-time evolution. If the center of mass of the system is allowed to 



drift, Eq. (7.96) will then contain an unphysical contribution that is usually difficult 
to consider. 



7.1.6 The Worm algorithm 

Here we present the Worm algorithm in the Canonical ensemble. This method of- 
fers an enhanced permutation sampling and also a way to compute, within the same 
simulation, both diagonal and off-diagonal properties of the system under study. Dif- 
ferently from the worm in the Grand Canonical ensemble, this method can also be 
applied in PIGS without any further adaptation: the Worm algorithm in the Canon- 
ical ensemble has not moves that create, destroy or change the length in imaginary- 
time of the polymers; these moves, in fact, would not be of easy interpretation in the 
context of quantum evolution in imaginary-time; Canonical Worm is based on moves 
that in PIGS can be applied at the time-slice at position M/2, and in PIMC can be 
applied anywhere in the path integral. The space of configurations that is sampled 
by Metropolis is enlarged by including also configurations with one open polymer, 
see Fig. |7.13| this polymer is called "worm" . A configuration with a worm is called 
"off-diagonal" (in worm notation: G sector) whereas a configuration without worm is 
a "diagonal" (Z sector) configuration. Starting from a polymer i, defined by the set 

of beads {r^} -_ v where eventually r/ = rj^, a worm at position m is constructed 
with the following operations: 

• Remove the kinetic correlation between r/™ -1 and r™ 

• Add a bead that is linked by a kinetic correlation only with its previous 
bead, f." 1 " 1 

• The beads and r™ are the two worm extremities on which off-diagonal 
properties such as the one body density matrix can be computed. 
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Figure 7.13: The worm in PIMC (a) and in PIGS (b). The gray lines represent the 
removed kinetic correlations. The worm in PIGS originates two half polymers that 
are not connected each other; in PIMC, instead, it originates an open polymer. 



In PIMC the timeslice m can be anywhere between 1 and M because any time sector 
is an equivalent representation of the system; in PIGS m should be in the time sector 
range in which the trial wave function can be considered an accurate representation 
of the ground state, in our implementation of the worm algorithm in PIGS, m is fixed 
at the central timeslice. 

The implemented worm algorithm consists of two Metropolis moves plus an exten- 
sion of the Brownian bridge and an extension of the translation move which deal with 
the presence of the worm extremities. There are two input parameters: C and s. The 
parameter C sets the ratio g/z between the number of G and Z configurations that 
are sampled during the simulation, g + z are the total MCS after the equilibration. 
There is not an universal relation between C and g/z but large values of C result in 
simulations with more off-diagonal sampling. Usually C is of the order of unity, but 
the best choice may vary drastically with the system under study. The parameter s 
is a tuning for the worm moves and specifies how many time-slices are to be involved 
by the swap and the open/close moves. 

Open/Close These moves allow to switch from the Z sector to the G sector and 
vice-versa. The Open move creates a worm from a diagonal configuration while the 
Close move closes a worm and gives a diagonal configuration as a result. In order to 
maintain the detailed balance of the sampling, these moves are coupled meaning that 
in a MCS there is always one attempt to Open/Close and whether to Open or to Close 
is decided with a random number: the Open and the Close moves should always be 
equally probable, so, for instance, if a random number uniformly distributed between 
and 1 is lesser than 0.5 in that MCS an Open move will be tried, otherwise a Close 
will be attempted, no matter whether the system is in G sector or in Z sector. The 
Open move is as follows 
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Figure 7.14: Panel a) The "open" move in the worm algorithm. The gray area on the 
left side represents the Gaussian probability distribution from which the coordinates 
of t\ are sampled. Panel b) The "close" move in the worm algorithm. Grey beads 
and lines represent the removed segment of the polymer. 



If there is already an open polymer, reject the move. Otherwise select a random 
integer number % between 1 and N, a random integer number j between 1 and 
s and a random integer number m between 1 and M. 

Create a worm in the i-th polymer at the m-th bead, the added bead is 
sampled from the probability distribution of a free particle propagator with 
time-step jSr 



Go {fr j ,r*j8r) 



(A7T\j5ry d/2 



iXjSr 



(7.97) 



• From the bead r™ J remove the beads r™ 3+ ,r^ m J+2 ,...,r™' 1 and build a 
discrete free-particle path to the newly created bead r". 

The probability to accept this move is 
p = min{l,P } 

nc E7=m-j [u [RT W , gg) - u {Rf d - RfA)} 

f = j ; r . 

The Close move is as follows 

• If the configuration is diagonal, reject the move. Otherwise there is a worm, 
say in the polymer % at the bead m. Select a random integer number j between 
1 and s 

• Remove the worm extremity that is linked by a kinetic term only to its 
previous bead r/" -1 
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• Replace the beads between r™ 3 and r™, with a free particle path 
The probability to accept the move is 



V m 

Pc = ~Wr E P ( R i ew ^ ~ u ( R ° ld - R t d i)] 



p c = min{l,P c } 
V 

NC 

l=m-j 

Gofrr^r^jSr) . (7.99) 

These moves can be optimized in term of computing efficiency with an automatic 
rejection whenever the distance between r/™ - -' and f ^ is such that the quantity 



Pi 

2 



AXjSr 



(7.100) 



is larger than some arbitrary quantity of order unity. Our choice was to set it equal to 
4. This avoid very small acceptation rates when the worm extremities are far away. 
In order to maintain the detailed balance, this kinetic test must be applied both to 
the Open and the Close moves. 



Swap This move is attempted only in the off-diagonal sector and implements the 
sampling of permutations. Consider a configuration in G sector with a worm in the 
i-th polymer at bead m. 

• Select a random integer j between 1 and s 

• Select a bead r™ +3+l with probability 

N 

£T = ^Go(r7,C +J+1 ,j'H . 

n=l 



(7.101) 

(7.102) 



• Evaluate the quantity 

N 

Z K = Y,G {?r + \?n +3+ \3?>T) • (7.103) 
n=l 



• Consider the bead r* m 3 and insert a new bead r*X = r J™ that is connected 

»* Pi k Pi k 

by a kinetic term only to its previous bead. 

• Replace j beads of the polymer i k , namely r™ +1 , r™ +2 , r™ +3 with a Brow- 
nian bridge starting from r" and ending at r™ +3+1 . With this operation, the 
bead f t u is no longer a worm end: the Brownian bridge swaps the polymer % 
with the polymer % k and the new worm extremities become and r™. 
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The probability to accept a swap move is 



mm 



i^ e i:7Jz[u{Rr w ,R?+r)M R f d ' R t&)}\ . (7.104) 



A remark is necessary here: the interpolymer correlations in the time sector (m, m+1), 
namely U(R m , R m +i) must take into account the worm extremities fp. and r™. The 
worm extremities contribute to the interpolymer correlations like the other beads but 
with a weight factor of 0.5. Using a symmetrized form for the density matrix, such 



as Eq. (B.l) automatically gives the correct weights: 



U(R m ,R m+1 )p A = e~^^<k< r hk) e ilsrT.h{^ e -iFT. h<k v(r™ +1 ) 
U(R m -i,R m )p A = e'^^^^^^e^^^' 1 ' 7 ' 1 ^ e ~^^ h < kV ^ 



(7.105) 



where 



= ifM ^ (7.106) 

^ if // Pi K J 



h i 



and 

* = «MJ . ( , 107 ) 

I \ r h ~ r fc I > it h = Pi 

The swap moves are very efficient in the permutation sampling for two reasons: first, a 
permutation that involves several polymers is obtained with a certain number of swaps 
that are more likely to be accepted; second, the Worm itself, being an open polymer, 
has a better probability to avoid overlaps during swap moves. The mechanism used 



by the Worm algorithm to sample the permutations space is depicted in Fig. 7.15 
the basic idea is that two topologically different diagonal configurations are connected 
by at least three successful Worm moves: the open move generates an off-diagonal 
configuration; the swap move (or a series of swap moves) samples the permutations 
space and, finally, the close move returns the system to a diagonal configuration. 

Brownian bridge extension The Brownian bridge in a worm computation re- 
mains unchanged; however, if the attempted reconstruction involves the worm ends, 
i.e. starts from a bead ji that is before the imaginary-time position m of the worm 
and ends to a bead j 2 > m, the reconstruction of the segment is split and the position 
of the worm extremities is updated: 

• The first worm extremity, is updated from the distribution probability 



G (r?\r*, (to - ji)5r) 



(47rA(m-j 1 )5r)^ /2 
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4A(m. — ji )5t 



(7.108) 




Figure 7.15: The sampling of permutations in the Worm algorithm: from an off- 
diagonal configuration (a), the first swap generates an open polymer of length 2/3 (b). 
Another successful swap results in a configuration with a polymer of length 3/3 (c); at 
this point, a successful close move yield a diagonal configuration with a ring polymer 
of length 3/3 (d). Grey beads and lines represent the removed segment of the polymer. 



The beads rr ' ~,r 



31+1 -*j!+2 
i i i i 



T 



m—1 



are updated with a Brownian bridge that starts 



from r? 1 and ends to the freshly updated bead 

The second worm extremity, n m is updated from the distribution probability 



G (f*,rf ,(j 2 -m)<5r) 



(47rA(j 2 - m)5r)~ d/2 



Of-**) 

4\(j 2 -m)l 



(7.109) 



tm+l z?m+2 ffh-1 

from the second worm extremity r™ and ends at f" 



The beads fr + \fr +A fP~" are updated with a Brownian bridge that starts 



The new position of every bead is sampled here with a free particle propagator, as 
consequence the probability to accept this move is the same of that of a Brown- 
ian bridge between ji and j 2 that is expressed in Eq. (7.57) with U(R m _i,m) and 



U(R m , Rm+i) defined as in Eq. (7.105). 



Translation move extension The translation move for polymers with a worm 
has been extended as follows. In the PIGS polymer i with Worm extremities 

r} ,...,r," , r™ , r* M ). Here, 
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J and r, m is defined by a set of coordinates S 

l old l old •> 



for a polymer without worm, we define r v = f m . The translation of this polymer 
is parameterized by two displacement vectors di,d 2 ; from S, the move will generate 
a new polymer S new = fry 1 = rA, + d x , , + d\,rP). + d 2 , ■~,r i A f + d 2 ). The 

A ° v Lnew ''old 1 1 & old ' l old ' 1 & old ' 

probability to accept the move becomes 

a({R} ne J=mm(l,P tr ) 

Ptr = A l -A 2 (7.110) 



where A\ and A 2 are the probabilities to accept the translations of the two half 
polymers S x = (fg H , ...,r£J and S 2 = {r™ d , ...,f^); namely: 



k^i 

* T (fl loI<l )e-^ v ( Bl ^e-* TV ( Jia ow)...e-T :v ( J4 'o I d) ' 

p-^ ! V r (iZ T n neu ,) f> -(5rV(ii( m+1 ) ) p-y^Mn™)^ fD 

J _ 1 \ ±ll M new ) ,„ 113) 

2 " e-^ y ^o id ) e - 5ry(i? ( m + 1 ) oH ) ...e-^ y ^o id )^ r ( J R Mo! J ' 1 ' 1 

In the PIMC case, a ring polymer % with Worm extremities r 4 u and ry m becomes 
an open polymer. The translation in this case has only one parameter d, so that 
if the polymer is defined by S = (r^ 1 , , r™ , ...,f^, the new polymer will be 
S new = fry 1 + ri, n 1 " +d, r* m + ri, f- M + a) The probability to accept this move 

v l old 1 1 % old 1 ^old 1 1 lold 71 J r 

is 

a({ J R} ne J=min(l,P tr ) 

e -5rV'(iJi nero )_^ e -^V'(i?„ nem ) e -^V(iJ mnero ) e -(5rV r (iJ( m+1 ), lelu )^^ e -(5Ty(_R A / nem ) 



tr 



(.7.114) 



Particular care must be taken when applying the translation move to a PIMC 
configuration that has permutations. In this case, all the polymers that contribute to 
a permutation loop are translated by the same displacement vector d; the probability 



to accept such a move is of the form of Eq. (7.114) but the correlation V(R m ) must 
take into account not only the z-th polymer but also the other polymers involved 
in the translation: let the polymers in the permutation loop be (z'i, i 2 , ...in) and 
the remaining polymers the elements of the set W rem ; then, considering that the 
translation does not change the correlations between these polymers, 

H 

V{Rm):= £ E^r-^TI) • (7.H5) 

fc6W rem 1 = 1 

Normalization of the One Body Density Matrix The Worm algorithm pro- 
vides also a way to compute the correct normalization of the OBDM. We focus here 
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on the PIMC case, the PIGS case is analogous; the OBDM is defined by Eq. (7.89). 



This probability density is sampled in the G sector. The normalization of Eq. (7.89) 
is the partition function Z; Z corresponds to Eq. (7.89) with r = f 
homogeneous system, this is equivalent to the sampling of p\(r = \ f 



in the case of an 
■r | =0); within 



the Worm algorithm this quantity is related to the probability to switch from Z to 
G. For an homogeneous system, the OBDM pi(r) is the histogram of the distance 
between the two worm extremities, r? and fl TO , normalized as follows: 



(S(r-\rr-fT\)) 
V she ii(r)zCp 



(7.116) 



where C is the Worm parameter previously introduced, p is the density of the system, 
z is the number of Monte Carlo steps in the Z sector. The Monte Carlo average (...) in 
this context is the histogram of the distance \r" — r^ m |, each bin of the histogram has a 
width dr and is divided by the volume of the spheric shell V s h e ii(r) = V{r + dr) — V(r). 



The introduction of a repulsive factor such as (7.90) does not preserve the Worm 
normalization. This happens because the repulsive factor interferes with the prob- 
ability to switch from the Z sector to the G sector and vice versa. A workaround 
is to use a repulsive factor f w (r) that goes to unity for r = 0: with this repulsive 
factor, in fact, when f = f one obtains again the partition function Z. A straight- 



forward adaptation of the repulsive factor (7.90) that has been used in this work is 
the following 



fu 



1 + A + C 



1 + Ae~ Br2 + Ce~ Dr 



(7.117) 
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Appendix A 
Estimators 



Here follows the derivation of some estimators with a fourth-order approximation of 
the propagator. 

We consider the Pair Suzuki (PS) approximation: 



G (R m , Rm+i, t) 



(AttXt) 



Nd 



N 



Y[dr*exp 



At £->i=l 



i=i 



exp 



■ | Ef<, [ve (r%) +4v c (rt j ) +v e (r^ 1 )] ^ 



where rV\ = Iff 1 — rTl is the distance between the i-th bead and the ?'-th bead 
at an imaginary-time defined by the index m; iV is the particles number, d is the 
dimensionality of the system, 2M is the effective^] beads number, A 
and 



2m' 



T 



J_ 
2M 



v e (r) = v (r) + -ar A 
-y c ( r ) = u (r) + i (1 - a) r 2 A 



9v (r) 

dr 
dv (r) 

dr 



(A.2) 
(A.3) 



We remark that (A.l) is the Green's function that involves two adjacent real times- 
liecs, thus the integration variables {r* } represent the fictitious bead required by the 
PS approximation. 



A. 0.7 Total Energy 

The thermodynamic estimator for the total energy is defined as follows 

e -4m% < a - 4 > 



1 With the parameter a set to zero, the odd timeslices are those which describe the system, the 
even timeslices are the fictious beads used to express the fourth-order approximation of the Green's 
function 
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where Z is the partition function, 

M-l 



Z = tr{p}= I dR 1 ...dR M Y[G(R m ,R m+1 ,r) 

J m=l 



(A.5) 



The imaginary-time derivative applied to the productory in (A.5) yelds a sum of M 



terms that may be viewed as the energy evaluated at an imaginary-time sector. The r 



derivative applied to ( A.l ) gives three terms: the first comes from the normalization of 
the kinetic part, the second from the gaussian which expresses the kinetic propagator 
and the last one from the term involving the inter-polymer correlations. After some 
straightforward algebra, one finds 

Nd 1 ^(rf -r7) 2 +(r*-rf +1 ) 2 



L'l" _ I _ _ \ 

^ ~ ' 2r 4Ar 2 ^ 2 

i=l 

A. 0.8 Kinetic Energy 

The thermodynamic estimator for the kinetic energy is defined by the following for- 
mula 

l3Zdm fiZd\ K ' } 

The arguments of the previous paragraph apply here too, resulting in the following 
expression 

Rm _/Nd 1 " (rf - rtf + (g - rT +1 ) 2 
\ 2r 4Ar 2 ^2 



i=l 



(A.8) 



A. 0.9 Pressure 



The thermodynamic estimator for the pressure is obtained from a volume derivative 
of the partition function: 
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In order to compute this volume derivative from ( A. 5 ) one has to perform the following 
change of variables 



dr = dxdydz = Vdf 



(A.10) 



The estimator is composed of three terms, the first arises from a factor V 2NM coming 
from the Jacobian transformation of the differentials dr? 1 , the second is the derivative 
applied to the kinetic factor of the propagator and the last one comes from the inter- 
polymer correlations part. Shifting back to the former integration variables, one 
obtains 



pr< 




S ^( rT _ r7 )2 + ( r7 _ rT+ l) 

^ 2 



6Vd ^ ■ 13 



Or 



+ 



dr 



+ 



dve 
dr 



,m+l 



(A.ll) 



A.0.10 T=0 limit 

Even though the previously introduced estimators are derived from a finite- 
temperature background, it can be shown that they are valid also in the zero 
temperature limit. Let's show this for the Hamiltonian operator 

^oH^o) = - hm 4 log ! dR 1 dB? M ^T (R 1 ) G (R\R 2M ,(3) ^ t (R 2M ) = 

I /3^oo dp J / v / 



lim — r / dR l dR 2M ^! T (R 1 ) 



dG (R\R 2M ,p) 



dp 



^ T (R 2M }A.12) 



Because of the formal similarities between PIGS and PIMC, chosen a large enough 
imaginary-time (3, this expression evaluated at the central timeslices expresses a zero- 
temperature quantum average. 



A. 0.11 Virial Energy Estimator 



Eq. (A. 6), such as any estimator involving the Kinetic Energy, contains an higly 



fluctuating term which comes from the derivative applied to the kinetic part of the 
propagator. This results in a variance of the averages that increases with the beads 
number M. The virial estimator provides a way to get rid of these fluctuations. Let's 
derive it for the total energy estimator. Consider the quantity 



E 1 



L+l 



I NLd 




(A.13) 
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where 



a 



m=l 
L N 



A\T 2 M 



& = E E a [«• + 4«o (4) + 



(A.14) 
(A.15) 



m=l i<j 



The quantity represents L times the total energy of the system, where L is 

a parameter arbitrarily chosen between 1 and M. For simplicity, let's rewrite the 
definitions of a and U in a more treatable way: 



m=l 



4Ar 2 M 



2L 



£7 = ?7 (iT\ R m+1 ) 



(A.16) 
(A.17) 



m=l 



where now the index m denotes every beads, both physical and fictitious, and 

U* (r m 1 r m+1 ) = v e (r™) + 2v c (r™ +1 ) m odd (A. 18) 



N 



U (R m , R m+1 ) =J2U*{ 



i<j 



Now define dR m = ]jf =1 dff, (R m - R n ) = £f =1 (ff - fj 1 ), ., 
consider the quantity 



(A.20) 

EJV 



* a and 



G 



SdR\..dR^Y^ = 2(R m -R l )(- 



1 ' 9 PYn"^ 

-ft? 



(A.21) 
(A.22) 

m Dm— 1 



with g = a+^-. If we make a change of integration variables, namely S m = R m —R 
eq. (A.21 ) becomes 



G 



JdS\..d6* L E 2 m L =2 6 m (- 



I ggm e^jj 



/diP...di^£^ 2 exp 



-ft? 



(A.23) 
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vanishes if r\ << V<*, thus 



The integral at the numerator of (A. 23) can be computed by parts. The surface term 

N(2L-l)d 



G 



n 2L 

IE 



m=2 



dS m 
85^ 







(A.24) 



The quantity in the RHS of (A. 21) can be expressed by explicitely computing the 
derivative over the positions 



2L 



g = ( i Rm - Rl ) 



\m=2 



da 
dR™ 



vh 1 )dRm l 



(A.25) 



The first term on the RHS, after some algebra, becomes 

21 



da 



E (*■ - Rl ) ^ = 2a + ivsu ( fl2i - fl2i+1 ) ( R2M - R ') < A - 26 > 



m=2 



dR r - 



A\t 2 M 



Equating both (A.24) and (A.25) and using (A. 26), the quantity a may be re- 
expressed as 



a 



N(2L-l)d 



8Xt 2 M 



(R 2L - R 2L+1 ) (R 2L+1 - R 1 ) 



1 2L 



m=2 



9Z7 



(A.27) 



Substituting a in eq. (A. 13) we finally obtain the virial estimator for the total energy 
per particle: 



2L N 



d 



2t 4Xt 2 N 



1 N 
Tin ^ ( 



)(' 



2riV 



EE 



i=i 

0E7* (rf, r™ +1 ) 2 y> 0/ // . // 



2L 



m=2 i<jl 



dr. 



m=l 



dr 



(A.28) 



This estimator may be used also in the zero temperature limit if one performs the 

following substitutions: R 1 ^ R T , L ^ L ,Em=i ^ Y?m=r and Em=2 ^ Emt^+i, 
where T represents the index of the first time-sector that can be considered a ground- 
state description of the system and L is an arbitrary number between 1 and the 
number of physical timeslices available for the evaluation of ground-state expectation 
values. 
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Appendix B 



Higher order approximations for 
the density matrix 



In this appendix we show approximations for the small imaginary-time density matrix 



that go beyond the Primitive Approximation (PA) introduced in Eq. (2.16). We have 



already shown the "Pair" Suzuki (PS) approximation in Sec. 2.1, we will show now 



the Pair Product approximation (PPA) and the Multi Product Expansion (MPE). 
The MPE and the PS, as well as the Primitive Approximation, have the advantage to 
be analytic and thus estimators can be derived exactly; the PPA on the other hand 
is numerical and only a restricted set of estimators, such as those diagonal on the 
coordinate representation and the one body density matrix, can be simply derived. 
The PPA, however, requires fewer imaginary-time projection than the other two and 
this feature could be useful in some contexts. 

As mentioned before, given the Hamiltonian H = T + V, one has to use a small 
imaginary-time approximation of the propagator e~ SrH in order to obtain an analytic 



expression for G (R, R', St) = \ R 



R'). The simplest approximation is the PA 



G 2 (R l} R J ,5r) = t r^e-^t r? v > 



(b.i; 



which is correct up to second-order in St; this approximation is obtained by ignoring 



the commutator 



T, V 



when factorizing the Hamiltonian. The effective potential 



U (Rm, Rm+i, St) for the beads represented with the PA is 

U (R m , R m+U St) = y [V (Rm) + V (R™+i)\ 
so that the density matrix takes the form 

G (R m , R m+1 ,5r) = G (R m , R m+ i, St) e -^.Vi,fr) 
where Gq is the density matrix for free particles 



(B.2) 



(B.3) 



Go {R m , Rm+i, St) = ( R 



-StT 



R 



m+l 



(B.4) 
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B.0.12 The Pair Product 



Following Ref. pQ, the PPA is a decomposition of the density matrix in which the 
effective potential is written as: 



St) = u 2 (?' 



^*m' '"m+l 



r 



3 

m+l> 



St) 



(B.5) 



i<j 



where is the position of the z-th particle at a timestep r m = mSr. u 2 is the exact 
effective potential for two atoms. This approximation states that, if the imaginary 
time is sufficiently small, the many-body propagator can be described as a product 
of two-body propagators. The density matrix for two particles can be obtained with 
different methods. For instance, one can use the matrix-squaring method; this is 
shown in detail in Ref. pQ and we limit here to state the final result for the two-body 
effective potential, 



u 2 {ri,r 2 ,ST) 



u Q (fi,Sr) + u (f 2 ,5r) 



(B.6) 



k=l j=0 



where q = (\r\ \ + |r*2|)/2, s = \r\ — r*2| and z = \r*i \ — \r*2\. The first term is the effective 
potential of the PA and the functions off-diagonal terms that are obtained 

from the partial wave expansion of the two-particles propagator. These off-diagonal 
terms are usually obtained by numerical means and an analytic description of the 
pair density matrix is not available. 



B.0.13 The Multi Product Expansion 



From Eq. (B.l), the 2n-th order multi-product expansion is built with the following 
relation 



G 2n (5t) = J2^G2 {Sr/kt 



i=i 



n 



h 2 — h 2 



(B.7) 
(B.8) 



Let's focus on the case n — 4. Following ref [2], the convenient choice for {fcj} 
that yelds an eight order multi-product approximation suitable for PIGS calculations 
is {ki} = {1,2,3,6}. This, infact, produces elements with time-steps Srjk\ 
that have a common divisor 5t/6; these elements can thus be represented by a path 



integral with a time-step St/6. This choice for {ki}, combined with Eq. B.7, gives 
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the definition of the propagator used in Sec. [5] 



G 8 ({Ri} 7 i=1 , Q5t) = G (1, 2, 5t) ...G (6, 7, 8r) 
~54 
35 

27 
~40 



x 



e -^Vi e ~ST% e ~STV 3e -STV4 e -8TV 5e -5rV6 e -^V7 



e -5rVi e -25rV 3 e -28rV 5 e -5rV 7 



2 

15 



_| e -^SrV le -3STV 4e -l5rV 7 



840 



_ e -35rVi e -35rVV 



where Gq has been defined in Eq. (B.4). 
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